{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assembly kmer estimates\n",
    "\n",
    "We wish to the following properties of a genome based on the kmer histogram. Formally, we have as our input:\n",
    "\n",
    "- $K_i$, where $i = 1..100001$, as the number of kmers at coverage $i$\n",
    "\n",
    "We want to infer:\n",
    "\n",
    "- $G_i$, where $i = 1..100$, as the number of kmers with copy number of $i$ in the genome sequence\n",
    "\n",
    "This $G_i$ spectrum is interesting since it allows us to estimate\n",
    "\n",
    "- Genome size $g = \\sum_{i=1}^{100}{i * G_i}$\n",
    "\n",
    "- Ploidy $p$ where $\\sum_{i=1}^{p}{G_i} > 0.95 * \\sum_{i=1}^{100}{G_i}$, i.e. with this ploidy you can explain most of the unique kmers in the genome\n",
    "\n",
    "- Repeat contents $R = \\sum_{i=p}^{100}{i * G_i} / \\sum_{i=1}^{100}{i * G_i}$, ratio of all high copy kmers (larger than ploidy $p$) over all genome kmers\n",
    "\n",
    "- Heterozygosity. Formula TBD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 342,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set_context(\"talk\")\n",
    "import matplotlib.pyplot as plt\n",
    "import logging\n",
    "logging.getLogger(\"matplotlib\").setLevel(logging.WARNING)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import nbinom\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nbinom_pmf(x, lambda_, rho):\n",
    "    n = lambda_ / (rho - 1)\n",
    "    p = 1 / rho\n",
    "    return nbinom.pmf(x, n, p)\n",
    "\n",
    "def plot_negative_binomial(x, lambda_, rho=1.2, color='b'):\n",
    "    y = nbinom_pmf(x, lambda_, rho)\n",
    "    plt.plot(x, y, '-', color=color)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 292,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbd19d2c090>]"
      ]
     },
     "execution_count": 292,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd3hUVfrHP2fSE0IKodfQsSCIgg3FhmXt2AvqWnZXXXvZddefrq666rqWdXV1de1rb9hXXGwgIEgR6RAIEEhI73XO7493LncIM8lkMuVOcj7PM8+9uXPPnTMJ3O9961FaawwGg8Fg6CiuaE/AYDAYDLGJERCDwWAwBIUREIPBYDAEhREQg8FgMASFERCDwWAwBEV8tCcQTpRSzYhIVkZ7LgaDwRBD9ATcWus2NUJ15TRepZQbUBkZGdGeisFgMMQMFRUVAFpr3aaXqktbIEBlRkZGRnl5ebTnYTAYDDFDZmYmFRUV7XpuTAzEYDAYDEFhBMRgMBgMQWEExGAwGAxBYQTEYDAYDEFhBMRgMBgMQWEExGAwGAxBYQSki9BY2EjJxyVU/VgV7akYDIZuQlevA+kWrH1oGwW3rgNAxykO/HESPcb3iPKsDAZDV8dYIDGOu8nNuv/bvOtn1aJZevXGKM7IYDB0F4yAxDirnikmrb6RFuCf8SMBaP6ulLKvyqI7MYPB0OUxAhLjrPpLAQBLk3tx7FMDWUFPANbeuzWa0zIYDN2AgAREKZWklHpAKVWglKpTSs1XSh0d4NiBSqk3lVLlSqlKpdT7SqlcH+dpP69fd/RLdRdKVtSRs1X6fKWeN5CZFyvmZg8AoPLrctxN7mhOz2AwdHECDaK/AMwAHgXWA5cAnyqljtBaf+9vkFKqBzAHSAfuBZqBG4CvlFITtNat/SyfA6+0OrYgwDl2O5a8XEE8UE0cZz+QRUICTLg8Cx6EhKYWKudXkjk1M9rTNBgMXZR2BUQpNRk4F7hBa/2o59hLwArgAeDwNoZfBYwEJmmtl3jGfuoZewPwf63OX621bi0gBj8UflXFQGBbj5706q0AOOyUJFY/mMYIatjybpkREIPBEDYCcWGdCTQBz1oHtNb1wHPAYUqp/u2MnW+Jh2fsauBL4GxfA5RSKUqp5ADm1e1xrZZuy82je+46duCBsCw+C4AdH5VGZV4Gg6F7EIiATEQsg+pWxxcCCpjga5BSygWMBxb5eHshMFopldrq+OVADVCnlFqulDo9gPl1S5pqWuhdKX+S3kek7zqemAh1e2cDEL++iqbSpqjMz2AwdH0CEZD+wHYfx61jA/yMywaS2hirPNe2mAfcDpwKXO0Z+65S6jx/E/ME5v2+gC67FOHKD6qJR1aTHH9Oz93eG3pSBs0oFFC5wKzmazAYwkMgApICNPg4Xu/1vr9xBDpWa32o1vpxrfWHWusngf2BTcCDSikVwDy7FetniTDscCUzZnLibu8ddkwcmxHjrvC71oajwWAwhIZABKQOsQZak+z1vr9xBDkWrXUN8E9gEDDGzzmZbb2ACn/Xj3WqFknPq7J+PWktr1OmwHolrUy2f2N6YxkMhvAQiIBsZ3dXk4V1rMDPuFLE+vA3VuPbveXNFs82u53zuh2J22tkO27PnlcpKVDRR+IiTSuNBWIwGMJDIAKyFBjrqenwZopnu8zXIK21G/gJOMDH21OAdVrr2nY+e7hnuzOAeXYb3C2a7Fox3rLG+/YgJu4lf67k0nqaK5ojNjeDwdB9CERA3gYSkAwpQCrTgUuBuVrrAs+xIUqpsT7GHqSUmug1dgxwFPCW17Gc1h+qlOqF1JHkaa3XBfyNugGFPzeQjFSZDzq4dSKb0PcQW++rlxorxGAwhJ52Cwm11guUUm8hwez+wAbgYmAoUpFu8RJwBJJdZfEkcAXwiVLqYaQS/UbEdfWI13nXKKVOBT4C8oGBwJVAH+C0oL5ZF2b9HLE+3MDoab4tkHEHxLOFFAZTR+WiKjKPMAWFBoMhtATaTHEm8Jhn+zhikZyotZ7b1iCtdRUwDfgOuAO4B3GJHaG1LvE6dR5QjIjNP4BrgSWe8z4K9Mt0FwoXi+evOC6ZzN6+/4T77gvr8QTSvzUWiMFgCD0B9cLyVJ7f4nn5O2ean+NbgbPauf5/gf8GMhcDVK+qIweo7OkvgxpycyE/oQc07aRqeU3kJmcwGLoNpp17DKK3iAXSMtB3/APA5QKGyvtqay3arSMxNYPB0I0wAhKDpJVJDCRllH8BAcgcL+/HNblp2OKrntNgMBiCxwhIjNFU56ZXowhI74n+XVgAgyan0OzJaahd3V7GtMFgMHQMIyAxxsbv6ojz7A+b2rYFMnKsiwJP0X/1ShMHMRgMocUISIyx9QdpI9aEYsTBvrrE2IwcCZtJA6D4R2OBGAyG0GIEJMYoXS2xjNK4JBKT2u4xOXw45HuaKlb+ZATEYDCEFiMgMUbtJrFAatPatj5AemJVZ4mANG80AmIwGEKLEZAYo7lALJCm7MAWbXTlioDEVzWZxaUMBkNIMQISY8SXiAUS1799CwQgYx870G4ysQwGQygxAhJjpFaLBZIyLDALZNje8ZSSAEDder/LrxgMBkOHMQISQ7hbNJnNIiBZYwKzQEaOhALPwo9GQAwGQygxAhJDbF/ZSIJnHfS++wZmgXgLSNnPRkAMBkPoMAISQ2xdbLcjGXpAYBbIiBHsKiasXFPfztkGg8EQOEZAYoiiFSIA1cSTPTigRsqkpUG1p2tvS76xQAwGQ+gwAhJDVK4TC6Q8KQnVdg3hbrgGiYDEVTXRXGWWtzUYDKHBCEgM0ZAvFkh9z8DiHxapI+2mi3UbjBViMBhCgxGQGEIXiQWiewcW/7DoOyaBWk8LxvoNJg5iMBhCgxGQGCKxQgQkYUDHBGRYrtoVSDcWiMFgCBVGQGKIlIZG2Q5I7NC4YcPsVN5aUwtiMBhChBGQGEFrTXqz9LJKH5LQobHeAlKxygiIwWAIDUZAYoTK7S0k4QYga3jHLJChQ72q0U0MxGAwhAgjIDHC9lWNu/Z7j+6YgKSmQk2GxED0jnrcje6Qzs1gMHRPjIDECCXr7Vbs/cZ1zIUFED9ULBCloX6zsUIMBkPnMQISI5RvFAukhjh6ZMe1c/aeZI5KohmpPjSZWAaDIRQEJCBKqSSl1ANKqQKlVJ1Sar5S6ugAxw5USr2plCpXSlUqpd5XSuW2M2aKUsqtlNJKqcxAPqerU71FBKQqoWPuK4uhw13s8KTymloQg8EQCgK1QF4AbgBeAa4D3MCnSqmD2xqklOoBzAGmAvcCdwL7A18ppbL8jFHA44BZ/ciL+u0iIPXJQQrIUEwtiMFgCCntCohSajJwLnCr1vpWrfUzwFFAPvBAO8OvAkYCJ2qtH9JaPwJMBwYiguSLiz1jngvsK3QPmnZKDKQ5vePxD4DBg70zsYyAGAyGzhOIBXIm0AQ8ax3QWtcjN/jDlFL92xk7X2u9xGvsauBL4OzWJyul0oH7gbuAsgDm1m1QZZ4srMzgLJBBg2wBqV5rBMRgMHSeQARkIrBaa13d6vhCQAETfA1SSrmA8cAiH28vBEYrpVJbHb8DqACeCmBe3Yr4ahGQuN7BCYi3BdKwqR6tdcjmZjAYuieBLCrRH9jm4/h2z3aAn3HZQJLXea3HKs+1NwAopUYh8ZUZWutmFUC/cqVUeTunZLR7kRghuU5cWMn9g3Nh5eRAcWIyNAL1bhq3N5LUwZ5aBoPB4E0gFkgK0ODjeL3X+/7G0YGxjwDfaK0/CmBO3Y70JrFAegwOzgJRCuIGmbbuBoMhdARigdQhlkRrkr3e9zeOQMYqpY4HjkfcZQGjtW4zxddjocS8FVJb1kIqLQBk5gYnIAB9h8ZRvDGRHBqp21BH5lSTIW0wGIInEAtkO+Jqao11rMDPuFLE+vA3VmO7tx4EZgFVSqlhSqlhgHV3G9JOoL7Ls32VXYWeMyo4FxbsHgcxtSAGg6GzBCIgS4GxnpoOb6Z4tst8DdJau4GfgAN8vD0FWKe1tmo9hgCnA3ler+u8rv9eAPPsshSvs/tg9R0bvAUimVieWpCNxoVlMBg6RyAurLeBm4HLgUdBKtOBS4G5WusCz7EhQKonTdd77P1KqYlWKq9SagxSR/IXr/MuAFo/Wp8LnANcCGzp4PfqUpRtaCQZqMdF5oCOtzGxGDwY5jmsFqRmZQ2b7t5E3fo6XAku+pzfh/6X9ycuJfjvaTAYIkO7AqK1XqCUegt40ONK2oAU+w0FLvE69SXgCMA7fepJ4ArgE6XUw0AzcCPiunrE6zM+bv25SikrPfhjrXV72VZdmqp8EZCquERcrvaz0/zhNBfWlke2sOGWDXjCOwBUzq9k69+2Mv7z8aSObp3lbTAYnESgrUxmAo95to8j1sKJWuu5bQ3SWlcB04DvkBqPexCX2BFa65Ig59ztqNsuMZC65ODjH2AJiLiwmoqbaK5s7vTcgqXw1UI23CjikTIyhdz7cul/ZX9UoqJ+Uz1Lpi6hekXr0iODweAkAnFhWZXnt3he/s6Z5uf4VuCsjk5Ma30XUpHe7WkqkhhIU4/g4x+wuwUC4sZKn5jeqWsGQ/XyalZfKp7OzCMzGf/peFxJ8izT/7L+LD9hOU1FTaw4ZQWTfpxEQmbnhNNgMIQH0849FigVAXFndE5AMjOhOTWBGiS+UL8x8m4srTXrrl2HbtKkjEph73f33iUeAD0n92S/2fvhSnZRn1fPmkvXmKp5g8GhGAGJAeIqxYUVl9O5J3GlYPAQFdWmijvf2UnF1xUAjPr7KJ/WRfrEdEb9YxQAxe8XU/SfoojO0WAwBIYRkBggqU4skKR+nbNAQFJ5t0eprbt2a/L+kAdAr5N7kX1ctt9z+/+yPzkzcgBYf+N6mkqb/J5rMBiigxGQGCCtUQQkbVDnBSSabd1LPimhztMJeNCfcrnpJujVC9LT4fjjYUurZO1Rj40iLj2OpqIm8v6YF9G5GgyG9jEC4nAa69z01JItlTGs88HkaKbybn10KwDpR2bxi6t78Le/QWkpVFfD55/DhAmwYIF9ftLAJIb9aRgABc8UULvOrDFmMDgJIyAOp3CN7brpNTI0Fsg2a2nb/HrcTe5OXzMQqldUU/6llPN8nDyI77+XmMzvfw/PPiuWSGkpnHceVFXZ4wZeNZDkYcnQgrFCDAaHYQTE4RStsduY9OlEGxOLwYNhu5XK64b6zZGxQnY8v0N2BqZw16cS+3jwQbjvPrjsMvj6a0hKgrw8uPlme5wrycWwe4YBsPPNnVQtrsJgMDgDIyAOp3SDZylbFL2HB1S20yaDB0MhSTR7GgZEwo2lW/SuTKqPW/qiURx4INzgtajx3nvDvffK/jPPwPLl9nt9z+9L2vg0ADb+bmPY52swGALDCIjDqdwsFkilK4G4uODbmFgMGgRuXBR6uuxHIpBe9mUZjTvke7yyoy8ATzwBca3aXV1/PYwdK/sPPmgfVy7F8PuHy7Vml1E6uzTsczYYDO1jBMTh1G6TG29tUufdVwA9e8prVyZWBLryFr5SCMC27Ax2kMIhh8DkyXueFxcHt90m+6+/Lu4si+wTssmYKku75P0+zxQXGgwOwAiIw2kqFAFpSAuNgEBkM7HcDW6K3y8G4M3SPoBYGv44/3yxklpa4LHH7ONKKYb/RayQqkVVlMyKbis17dZULqxk6xNbWXf9Otb8ag3rb17P1r9vpfKHStzNkUlOMBiiSeed6oaw4vYU0Ll7hq4f1ODBUPBzZIoJy+aU0VLVglbwre7NwIFw+un+z09MhGuugd/9Dl59FR56CBI8Xz3jkAyyT8im9NNS8u7Io9fJvVCd6E4cDM1VzWx7YhsFTxfQsNnXas1CQk4Cvc/uzcDfDiRtbFoEZ2gwRA5jgTgcV4VYIK5e4bFA6jbWhdUdZFkfG1MyKCORc86B+HYeWy68UFJ8i4vhs892fy/3nlwAan6qoejNyLU40Vqz46UdLBi+gLzb83aJR9LQJLJPzKb3Ob3JPj5bUo6RbscFTxbww14/8PPZP1O11GSPGboexgJxOEk1IiCJfUMrIF96BMRd46axsJGkfr6Wru8c2q0p+UBcTZ/XSluSswLoyzxwIBx9NMyeDS+9BCefbL+XPimdnDNyKH63mE13bqL3mb1xxYf3OaipvInVF6/e5TZzpbgYcNUA+l/Wn7Rxe1oX9ZvrKXqriIKnCqjfWM/Ot3ay862d9Dm3D8MfGk7yoOSwztdgiBTGAnE4qY3iwkodGDoXlnc/LAhfV97KBZW7sq/m0ovBg2HKlHYGeZg5U7azZkF5q+XEcu/OBQV1a+sofLkwhDPek5pVNSw+YPEu8ciZkcOUDVMY+deRPsUDIHloMkNuHsLkNZMZ98o4UsfJwlhFrxexcMxCNt+3mZb6Fp9jDYZYwgiIg2lp1vR0yw2459DQWiB1xFPmWUW4bn144iCln0i67Y7kVApI5cwzxTUVCKefDsnJ0NgIn3yy+3tpe6fR53wJyG+6cxMtdeG5GVfMrWDJoUuo31CPK9nFmOfHsM/b+5DUPzBrzRXvou8FfTlwxYGMeX4MCX0ScNe6yftDHj/s8wPFs4pNNpkhpjEC4mCKNzZhlUpkh6CNicXgwbLd6nFj1a4NT4+p0v+KgHxdL5Xnp54a+NgePeCYY2T/gw/2fD/3T7moBEXDloZdPbZCSdmcMpYdu4zmsmYS+iYw8buJ9L+kf1DXUi5F/0v6M2XtFAbdOAgVr6jfUM+KU1ew/Pjl1KysCfHsDYbIYATEwRSutvtg9RkTOgEZNEi2WxDXitUhN5Q0lTZR9YMEjheTTY8ecPDBHbuGJTiffiqWiDcpI1IYeO1AAPLvy6dhh/+MqI5SNqeMn37xE+46N8nDk9l/3v6kT+r8yo3xGfGMfHgkByw/gKzpWfJZ/y3jh/E/sO7adTQWNrZzBYPBWRgBcTAl6zwrEQL9Rocu3yEtDbKybAGpXRN6C6TsyzLQ0OxSLCeDI4+UFN2OcNJJ4vKqqoKvvtrz/aF/HEp8r3haqlvYcNOG0Mz7f17iMSKZCV9NIGV4SvsDO0DauDTGfzaefWbtQ8rIFGiBbX/fxvzc+Wy4dQONO42QGGIDIyAOpmKT3EiqVAIJyaH9Uw0eDFusVN51dWh3aH3xZV+UAbAyPpMG4pg+vePX6NfPDrr7cmMlZCYw4sERABT9p2iXyyxYyv5Xxk8n7S4eyYPDkzGllCLn5BwOXHEgwx8aTnyveNx1brY8tIX5w+az5so1VC+rDstnGwyhwgiIg6nZKgJS09FH9wAYNMi2QNx1bhq2hs4FpLXedTOf1yiummAEBOwU3s8/9/1+v0v7kXG4tDhZc8UamsqCW7mw9IvSPcUjAum2riQXQ24ewkF5B5F7by7xWfG4a91s/9d2Fk1YxI+H/Mi2f2wLqYvOYAgVRkAcTGOR3AzrU0OXwmthFRNaDTdC6caqW1e3q9BuEdkMGQKjRgV3rWOPle2GDbv3xrJQSjH66dG4Ulw05Dew5oo1Hc5sKnq7aA+3VaRrNeLT4xl6+1AO2nQQI/8+clfqb+X3lay7Zh3fD/yepUcuZevjWyPWgt9gaA8jIA6mpVgskJb08FggzbgoTZIbZSgFxLI+apIS2EgaRxwRePpua/bfX+I1AF984fuctLFpjHx8JADF7xSTf39+wNcveKaAlWevRDdpUvdKZeLXE6Na6BffM55B1wziwJ8PZL8v96PfJf2Iy4gDN5R/Vc7669Yzf9h8fpjwA3l35lH1Y5VJBTZEDSMgDkaVi4Co7NALiJXKuysTa03oMrHK/ivxj8Vko1Ecdljw14qLs9N5/QkIQP/L+tP3QmkVn/eHPLY9ta3N67qb3Wy8fSNrf7UWNPQ8qCcTv51I0sDQV+QHg1KKrKOyGPv8WA4tPJR9PtyHfpf2IyFHrNGaZTVsvnsziyctZv7Q+eTdkUf9FmOZGCJLQAKilEpSSj2glCpQStUppeYrpY4OcOxApdSbSqlypVSlUup9pVRuq3OylVIvKqVWKaWqlFIVSqkflFIXKRXss2vsk1gtLqyEvuFxYQGsb/BkYoWoFsTd6KZ8jpSOf9cgpsPUqZ27puXG+vJL6dLrC6UUY54bQ9ax8pnrrlrHumvX0VK754CqxVUsPXzpLksla3oW+83ej4Ts0P+eQ4EryUXOSTmM/fdYDtlxCBO+mcCgmwZJBhfQsKWBzX/ezPxh8/np1J+o+L4iyjM2dBcCzQ19AZgBPAqsBy4BPlVKHaG1/t7fIKVUD2AOkA7cCzQDNwBfKaUmaK3LPKf2BIYD7wH5QBxwDPASMAr4vw59qy5CaoNYICkDwuPCAjsTq3Z1aASkcn4lLdVy015MFr162YtEBYslIGVlsGQJHHCA7/NciS72fndvVp6zktJPStn2920UvVFEv5n9SB2XSlNpE6WfllL+P7s3yuBbBpN7by6uhNgwxlWcInNqJplTMxnx0AhqV9VS+J9Ctj+7nabCJkpmlVAyq4Rep/Zi5KMjSRkW2hRkg8GbdgVEKTUZOBe4QWv9qOfYS8AK4AHg8DaGXwWMBCZprZd4xn7qGXsDHmHQWm8CWj+n/kMpNQu4Xil1p+5mjl63W5Pe4mljMiR8ArIZ6efUkN9Ac1Uz8emdqzex4h+lWWmUliVx6mHBxz8shg2T16ZN8O23/gUEIL5HPPvO2pdN92wi/y/5NBU1seWvW/Y4L23fNEY+OpKso7I6N7koopQiba80hv95OMP+bxjF7xeT/2A+1YurKfmghLL/ljH8L8MZ+NuBdGND3hBGAnnsOhNoAp61Dmit64HngMOUUm31dzgTmG+Jh2fsauBL4OwAPnszkAY407cQRiq2tZCIaGZmbui/fmoqZGfDJuyGgLUrO2+FWPGPH9xyY+5M/MMbyw32zTftn6viFLl35XLQhoMY8rshZEzNIGlQEmn7ptHngj7s88E+HLDkgJgWj9a4El30ObsPk36YxF5v7EXigETcdW7WX7eeFaetoLmyOdpTNHRBAnncnAis1lq3rmpaCChgArC99SCllAsYDzzj45oLgWOVUqla61qvMclAD8/rcOBS4Duttc/SXKVUua/jXmS0875j2b7K/sqhbGPizeDBsKw0gcb0RBKrGqlZUUPPKT2Dvl5TaRNVi6R9yewK6X91yCEhmSpTp8LLL8N334HWgVk1SQOTdq2lHkry8mDBAtnff38YPTrkHxE0Sin6nN2H7OOzWffbdRS+VEjJrBKWHLaEfT/al+QhppW8IXQEYoH0x4dAeB0b4GdcNpDUxljlubY3lwM7gTzgRWA+cEEAc+xyWG1MAPqNC5+AAJRlihVSs6JzTf2s9iU6QdqXxMfDxImdnaVgWSDFxbB6dWiu2VHy86W9yogRcN558hozRrLE1qyJzpz8Ed8znnEvjmPMs2NQ8Yqan2pYMnUJdXnhXYHS0L0IREBSAF9lsPVe7/sbRwfHvg8cC5wHvOI5lupvYlrrzLZeQMymo5RvkgysauJJzQhPgNeKg2xLCpGAeNxXZYMzaSSO8eMhJUQx3DFjoHdv2f/229BcsyPMnw8HHggffywWUM+ekJkp7335JUyevGfbeSfQ/7L+jP98PHE94mjIb2DpEUupzzfpvobQEMidqQ6xJFqT7PW+v3F0ZKzWeqvWerbW+nWt9UXAWmC2UqrbpZLUbBELpDohfOGfXam8zSIg1T8F33vJu33JsgSJLUye3Ln5eaNUx+IgoWT1ajjuOCgqkqLGN96AkhJ5ffgh9O8PlZWyhsl330V2boGQdVSWiEh6HA1bGlh+wvKgW74YDN4EIiDb2dPVhNexAj/jShHrw99YjW/3ljdvA4NpO9OrS1K/XQSkPiU87iuwBWRZlQhIU2FT0J1g69bW0ZAvxuZHOyT+EUoBAVtAImmBVFTAaaeJQAwYILGPs8+Wdd1dLnFpLVokqcqNjXLuxo2Rm1+gZBySwT4f7INKUNSurOXns35Gt3SrxEZDGAhEQJYCYz01Hd5Yi5Mu8zVIa+0GfgJ8JV1OAdZ5B9D9YFkeMRsMD5ZmTxuT5jC0MbGwXFhLSuxMrJqfg3NjWdaHKyeBHyvkeuESkPx8eUWC226T+EZiIrz7ru+eXgMGiGsrJ0eskksvBbd7z/OiTdaRUtkOUP5lOXl3+GguZjB0gEAE5G0kjfZy64BSKgnJkJqrtS7wHBuilGpdMvY2cJBSaqLX2DHAUcBbXsd6+/nsyxBL5ccA5tmlUJaLISv8Lqw64oj3ZOfULA9OQKz4R83YbECRnt75AsLW7LcfpHvWdYqEFbJgATzjySF84IG213MfPlyyxEBcbE8/Hf75BUPfC/oy6EZ5csi/P5+ST0qiPCNDLNOugGitFyA3+wc97UyuBP4HDAVu8zr1JWBVq+FPAhuBT5RSNyulrge+QFxXj3idd7VSaqlS6l6l1BVKqVuVUvOB04GntNbrg/2CsUp8tVgg8b3DZ4EMHGjvN+eKgVm1uKrD13E3uimbIwKyOl3iH/vvL32sQkl8vJ0WHO44iNsNv/mNBMwnTIBrrml/zPHHw8UXy/5tt0nMxIkMf2A4GVPFqF/9y9VmAStD0ASa3jMTeMyzfRyxSE7UWs9ta5DWugqYBnwH3AHcg7jEjtBaez/6fImk7s4EnvCc60YskAD+63Y9UurkP3VyGNqY7PqMFHG7AFT2l/qPqoUdF5DK+ZW4a8RnM6dCBGTSpNDMsTWRioO89560TQF48kkRr0D429+kQLOqCv785/DNrzO44l2Me3kccT3jaCpsYu2v10Z7SoYYJSAB0VrXa61v0Vr311ona60na61ntzpnmtZ6j/IuT2bVWVrrDK11utb6FK31xlbnfKu1Pl1rPVhrneQ57xCt9b+7WwsTix7N4sJKHxzeInzLjbXN4xuqXVPb4arl0s8k/pG2bxrf/CxJd+EWkFWrpCYkHLjdcPfdsn/yyR1byz07G/7wB9n/5z9lHRMnkjw0mVFPSECn+N1idr67M8ozMsQisdFBrptRXdJCKtKQMCM3fBYI2IH0dSpdSjt1x91YJR+LMek6OJsKT+XN/vuHcJJeTJ5sr60erpTZWbNg+XLZv6Y2ujwAACAASURBVPPOjo+/6ioYMgSamuCee0I7t1DS98K+ZJ8gGXPrrllHc4Vpd2LoGEZAHMj2lbZPuveo8AqIZYHkFcWTOkZqNqt+CFxA6vPrdwXetwzsBUBaWvArELZHcrJt3Xzvtw905/jb32R74onBWVLJyfDHP8r+q6/Clj17OToCpRSjnhyFK9VF4/ZGNt2zKdpTMsQYRkAcSPFau8ir316RsUC2bIH0A8WN1REBsbJ44rPiWVgjcZQJE0IfQPfGcimFQ0CWLbPjK9dfH/x1LroI+vWD5mZ45JH2z48WKcNSGPqHoQBs+/s26jaaVieGwDEC4kDKNooFUoeLjL5hvBPjtTKhl4BULqwMeHzJRyIg2cdns2iJ/HMKV/zDwhKQRYvETRRKnnhCtmPH2ishBkNysi1AzzwD5e21/Ywig24YRNLgJHSjZuPvHFgFaXAsRkAcSLXVxiQ+vNYHwFB5+KSoCJImigXRkN9A3ab2n0Rbalp2Lc6U/Yte/Oip1glX/MPCSuWtqxOLIVRUVIjLCSRtt7NLaPz61+LOq6mBF17o9PTCRlxKHLn3ySKhO9/aScW8mG0hZ4gwRkAcSF2BCEhtcuQEBKAkO534LMlXLfuizM8Ir/M/LsFd50YlKBr2y6bEk5g9YUI4ZmozYIAEqQHmzQvddd98U0QpNVVcUJ0lI8O+zpNPRr86vaZGqurXrpW2K970Pb8v6QeIBbr+xvV00+RHQwcxAuJAmnaKX6apR/jX0RowwI5X5G9VZB4lLWbLZrcvIEVvSKVc1vQsfs6XuSYkwLhx4ZmrN+GIg1hWwowZ0m03FFx9tWzXrZOuvZFGa/joIzjySGkEOXasdDbOyoJTT5UuwwDKpRjx8AgAqhZUsfNtk9ZraB8jIE6kTB4PdWb4LZD4eDsOsnkzZB8raZ1lX5ah3f6fQpurmin9ROo/+pzdZ5cradw4O802nIRaQNauta2ZSy4JzTUB9tkHDve0Ao10e5OiIpg+XWpZvvpq93hRba2kKx98sFhJ1dWQeXgmvU6RTLrNd29u8+9vMIAREEcSX+lpY5ITgTsxthtr0ybIOlYqyZtLmqle4r+9e8mHJbjr3ahERc6pObsEZPz4ME/WgxUH2bwZCvz1g+4AL70k2yFDYNq0zl/PmyuukO2HH0JpaWiv7Y+ffpK142d7yn2PO06q67dskW7BL71kx6peeUX6fG3eDMPuHAbI2jDFH4SpUtPQZTAC4kCS6+RRMalfZJaCtwRk82ZIGZ5Ccq40Vmyr0d7256QTf/YJ2cRnxO8qvNtvv7BOdRf77SeZTtB5K0RreP112b/wQmnTHkpOPx169JC4g/U54WTtWskg27JFgvhvvQWffSat5gcNgtxcsToWLYLHHxe348qV4uYqy0kn+xdihW6+Z7OJhRjaxAiIA0lrEgskbXBkLZDNm2Wbc4Y0yNrx4g6fboyan2t2ZV8N+NUA6ursJV0jJSCJifKEDZ0XkB9/tFuOnHNO567li7Q0OOss2X/xxdBf35viYjj2WHFfZWfD3Llw5pm+z1UKfvtb+O9/JXEgL0+Ep9f1wwCoXlK9q8uAweALIyAOo6HGTbqWlhIZwyIjIMOGyXbTJtn2/6WsAVa/oZ7yb/YsYNj2xDYAUkalkH1cNj//bGcYRcqFBbYbq7OZWJZVMHYs7Ltv567lj5kzZbtwYfjWdG9pgQsukLVS0tLg008DE/Rp02Q9k6QkCfZf/nBPMj2uTGOFGNrCCIjD2LHKjnTmhLmNiYVlgRQUSKA1ba80eh4kaUjbn9190ciG7Q3seHkHAAOvHohyqV3uq7595RUprED64sXQ0BDcNbSW9F0Q66OztR/+OPxw+/dsxVtCzYMPijUB8NxzHVvQa9o0e+2Tzz6Dr4fIZKsWVu1a68VgaI0REIdRtMZO0O87NrIxELcbtm6V/X6X9QOksKxmpb3I1Prr1+OucROfHU+/S+QcK4AeKfeVhSUgjY126/WOsnixvbrh2WeHZl6+cLlsK+Tll8VaCCWrVsFdd8n+VVcF54qbOVNcWgA3vZCJa5KkdG/5q0ObeRmijhEQh2G1MWlE0WtogItQdJLBg+0nbysO0ve8viTnJqMbNasuXIW70c2Ol3ew802pDxj5t5HEZ8j8Ip2BZdG3r6wECMHHQT74QLajR8Nee4VmXv6wBGTrVpgzJ3TXdbsl06uxUX4fDz4Y/LUeeEB+Dy0t8Hih5HeXzS6jern/jDxD98UIiMOo3CwCUhWXiMsVJn9KK5KSoL+EPXbFQeLS4hj70lhQEkyd13ceq2eK8z7zqEz6zhRfldZEPAPLG8sKCTYOMmuWbE85JTTzaYuRI+HQQ2U/lG6sl1+WYDmIGyotre3z2yIlRdJ64+Lg3a3Z1OZIh+atj24NwUwNXQ0jIA6jtkBiILVJkXFfWbTOxALIPCxzV4+k5nIJ7Gcdm8Ve/9kL5TFZtm6FMo+LPJoCEowFkpdni9+pp4ZuTm1x4YWy/eADqK/v/PVqauD222X/3HPh6KM7f82JE+GGG0CjeLZc2jUXvlpIw44gA02GLosREIfRVORxYaVFJoBu4UtAAIb+bigHbTqIkY+OZOzLYxn/+XgS+9pzs9xXCQnSIiPSWAKybVvH19348EPZ5uR0bNXBzjBjhsRDKivtgHdneOghSX5IThb3U6i4805xbX7S3Je6xHh0o6bgyRBUbBq6FEZAHIYu8bQxyYisgLRO5fUmeWgyg64bRL8L++2yPCysJ/i99opMC5PWjB8vNQzQcTeW5b466aTwrl/iTe/ecNRRsm9lfwVLSQk8/LDs33ij3WAyFPToAX/5CzQQx1uNAwEoeKqAlroQR/8NMY0REIcRVykuLFev6LuwAiFaGVgW8fF2umpH3Fjl5fD117IfifiHN1a21wcfSPffYPnrX6WHVUYG3HJLaObmzbnnSmfl9xlAs1I0FTdR+Eph6D/IELMYAXEYibVigST2i44La8uWjrUdj1YGljfBxEE+/VRWC0xKkoaDkeT008Xiqa6Gzz8P7ho7d0obEoCbboLMzNDNz8LlEiukjCRm6z6AFJGawkKDhREQh5HmWaghbWB0XFhNTbB9e5un7qK2ViqXIXoWCNgV6T/+GPgTvZW+e8wxnctaCoacHDvYHawb6/HH5feflQXXXRe6ubVm+nTpkfU+4saqWV5D5feBr1hp6NoYAXEQzY2adLe4sHoOjayAePvPfcVBfBGtFiatOegg2TY3S2FgezQ2igUCkcu+ao3lxpo1q+NurKoqe+nd3/42dGuX+EIpCc6voSerkQWntj25LXwfaIgpAhIQpVSSUuoBpVSBUqpOKTVfKRVQwqBSaqBS6k2lVLlSqlIp9b5SKrfVOYOVUncppRYqpcqUUsVKqTmBfkZXoWhdE1YsN3tEZGMgaWnyZAyBx0Es91W/ftCnT3jmFQg5OTBqlOwH4saaO1eyoEAC6NHgtNMkflNTY4tZoPzrXxLDSUmxK8fDyYEHSvbYLAYA0p2gcWdjO6MM3YFALZAXgBuAV4DrADfwqVKqzeRHpVQPYA4wFbgXuBPYH/hKKZXldeqpwK3AeuCPwD1AT2C2UioEi4vGBoWrvduYRD6lyXJjBSog0SwgbE1H4iCffSbbiRPtAspI06uXuM8A3ngj8HHNzfDYY7J/2WW26IebP/4R/kcfqpCU3h3/3hGZDzY4mnYFRCk1GTgXuFVrfavW+hngKCAfaC/z/CpgJHCi1vohrfUjwHRgICJIFnOAIVrr87XW/9BaPwYcAqxGxKRbULpR3FctQJ9RkbVAYPeFpQLBCQF0C+/OvO3FeC0BOf748M6pPawW7x9/HLgb66OPpHeXUuGNfbRmwgQ4+sQ4PkP6nxX8swDdYoLp3Z1ALJAzgSbgWeuA1roeeA44TCnV1jPcmcB8rfWuVnda69XAl8DZXsd+1lrvtvyZ1roB+AQYqpRKCWCeMU/lJk8bE1cC8YmRaWPiTUdSebWOfgqvN5YFUljYtgAWFNiWU7QF5NRTJRurpibwbKy//122J5wgrVEiye23226s+k31lH4WoeUVDY4lEAGZCKzWWrfuprYQUMAEX4OUUi5gPLDIx9sLgdFKqdR2PrsfUA2EoOmD86nZKgJSE42KPDomIFu2QEWF7DtBQPbeG9IlxtumG8uq/k5Pj1z1uT969bKLCt9+u/3zV66E//1P9q+5Jnzz8sehh0Lu1FQWId5nE0w3BCIg/QFfiZ3WsQF+xmUDSW2MVZ5r+0QpNRI4A3hb+0k89wTm/b6ADH/XdyIN20VAGiLcxsTCOwbSnhvIsj4SE6PTwqQ1cXGyrje0XZFuua+OOUbar0SbGTNk++GH7a9pYmVejRwpa5xHg9tvhw88/+VLPy2lLq8TlZCGmCcQAUkBfP3Trvd63984ghnrsUzeAmqA2wOYY5fAamPizoyOgFit0evq2q8FsQRkr72ccSOG9gPpLS22BRKtG3BrTjvN7o01e7b/8yoq7A6+V18d+nXbA+W446BuQi92kggaCp42/bG6M4H8M6xDLInWJHu9728cHR2rlIoDXgfGATO01n5vZVrrzLZeQIW/sU7EVS4C4sqJroCAvUa4P5wUQLewBGTZMokrtOaHH+zOwU4RkL59ZbVCaNuN9cIL8p3S0uCSSyIxM98oBTfd6uJDjxWy7V87aKk3/bG6K4EIyHZ8u5qsY/4eQUoR68PfWI1v99a/gF8AM7XWXwcwvy5Dco0ISFL/6AhIaqqd1hqogEzwGQGLDlZBYUuLiEVrrED12LG2u84JnHmmbN9/X4ocW6M1PPWU7F90UXjalnSEM8+ExX3704zCXdpE8bvF7Q8ydEkCEZClwFhPTYc3Ho8zy3wN0lq7gZ+AA3y8PQVYp7Wu9T6olHoIuBS4XmvdyV6lsUcPq43JEF9GW2SwMnvWr/d/Tk2N/b4TAugWWVkwbpzs+3JjOSV9tzWnny7b8nLfKxXOmwdr1sj+r38duXn5IyEBLrg2ie+QIpQt/zBurO5KIALyNpAAXG4dUEolITf6uVrrAs+xIUqpsT7GHqSUmug1dgxSR/KW94lKqVuAm4H7tNZ/D+K7xDT11W56aqkDyRweHQsEYMQI2bZlgfz0kx1kd5KAgP84SEkJLFwo+04TkAED7JUKfbmxnntOtpMmOef3feWV8HmCuLGq51VQ87MPn6Ghy9OugGitFyA3+wc97UyuBP4HDAVu8zr1JWBVq+FPAhuBT5RSNyulrge+QFxXj1gnKaVOBx4E1gGrlFIXtnpFuN1d5NmxyvZd9B7tbAGx3FeDBkkqqpPwFhDvTLLZs6VvV3KyHXNwEt5urOZm+3hVld1w8Ze/jPy8/JGTA/tenMkWTx7Mtn8aK6Q7Emgux0zgMc/2ccQiOVFrPbetQVrrKmAa8B1wB1JVvhQ4Qmtd4nWq9Vw1CnjZx6t3gPOMWQq9BGTAvtETkEBcWEuXytYpT8PeWBXpxcW7fwfLfTVtmvSQchpnnCHb4mL45hv7+JtvisswORnOPz86c/PHtdcpO5j+/A5aakwwvbsRkIBoreu11rdorftrrZO11pO11rNbnTNNa71H+bTWeqvW+iytdYbWOl1rfYrWemOrc+7SWqs2Xps69S1jgJJ1IiD1uMgYEKHl8XxgWSBlZXbGUmucVIHemrFj7SCz5cbS2g6gOyX7qjVDhtgLY3m7sSz31RlnRD943pp99oH6w/vSiELVtFD0RlG0p2SIMKadu0Oo2uxpYxKfiMsV+TYmFpaAgG83lttttwJxUgaWhctlFxRaAvLTT3Zdi9PiH95Ybqx335VMslWr7O9w2WXRm1db/OrWRL5CWjFveMS4sbobRkAcQt02EZDalOi5rwCys+0nXV9urI0b7RoLJ1ogsHtjRbDdV0OHOqNq3h9WVXphocz93/+Wn3NzxfXmRE44AZYMEjdW84oqqhZXRXlGhkhiBMQhNBeKgDSlR1dAoO04iBX/SE3d3VpxElYgfcUKCUJbAnLccVII51SGD4f995f9N96wK88vvTR6left4XLBL27tyUYkz2XjY8YK6U449J9lN6TUE0TPjr6AWIszrV2753veFehx0QvVtMmUKSIUbjd89RV8950cP+GEqE4rICwr5LXXoKhIvkc0K88D4ZJLFV8kixVS/EYhzRXN7YwwdBWMgDiEhCoRkIS+0ReQsZ5qHqt4zRsnB9AtevaUAC/A66/LOu/x8XbnWydjxUFKPZ3Sp0+HwYOjN59A6NEDBv+yL3W4iGt0s/2lwmhPyRAhjIA4hNR6EZCUgdEXECtOsGbNnl15Y0FAwHZjWdbHoYeGd+3wUDF6tC3g4NzgeWt+fVM8/6MvAGv+WoCfBtqGLoYREAfgdmsymqVpcfpQ5whIRYW4USzKymQ1PIgdAdm6VbZOzr5qzaBBsnW54OSTozuXQBk+HMoPFzdWfH4Nld9XRnlGhkhgBMQBFK1rIhF5Yuu9d/T6YFlYMRDY3Y1lWR9Kwb77RnZOHcUSELdbtrEiIFrb6dNutyQCxArn3ZnOamRVrxX3mWB6d8AIiAPYtsyuQh84PvoCkpZm+919CciIEfbqf05l9GjJFAPIyHC+xWSxYAHk5dk/B7JSoVM48khYPMCz5O2nRTSVNEV5RoZwYwTEAexcJe6rRhR9xzhjdSbLjbV6tX3MyS1MWqOUCCFAnz7OTt/15vnnZdtXwgm88077q0M6BaVgyu/7UE0ccW7N+n/siPaUDGHGCIgDqFgnAlIRnxTVKnRvfGViLV4s20mTIj+fjlJXZ7diqa1t+1ynUFcnWWMAM2fKdv16u/I/Fjj/l3F8k9QPgLzHTTC9q2MExAHUbhEBqUmNvvvKwjsTC+Qm/PPPsn+ArxVeHMY339hdbbdtg507ozufQHjvPVnaNj4ebr7ZLtSMJTdWair0vEBWJUstqaP4v+VRnpEhnBgBcQAtO0RAmjOdJyB5edDQIO4rKyAdCxaIVX1uua58LdTkNF54QbYnnyxuN6smJJYEBGDmnT34iQwAFt1pguldGSMgDsBVIgKiekc/hddir71k29IiFemLFsnPw4dLvyynYwlIbq5sZ8/2f64TyM+352hVnlsCsno1rFwZlWkFxZAhsGOyBNMTFxbT4HlAMnQ9jIA4gORqz1rog5xjgQwYYDdVXLHCFpBYcF9t2mQH/6303S+/jNp0AuLllyVY3qeP3XJl0iRpAAmxZ4Ucd28OFcQTpzU//MkE07sqRkAcQHqjp4gw1zkCopTdDiTWBMSyPnr1gosvlv2NG3dPj3USWtvuqwsvlDXHQf4GVm+st97yOdSxTD06jh/7SCyk7OUCdIsJpndFjIBEmariFtK1RHt7jXWOgIAtIEuX2k/0sSQg06fLU3xWlvzsVCtk7ly78/Gll+7+3tlny3bFitjKxlIKRt0gApJe08C6l0vaGWGIRYyARJmtS2z/cL99nCkgS5bYtQhWu3Gn0thoC8Xxx0vHYKuJolPjIFbtxwEH2L9zi8mT7fb6r7wS2Xl1ljOvT2Vpgqj3iru3Rnk2hnBgBCTKFP5sC8ig/ZwTRAf7Zmat5jd6tFR1O5m5c6G6WvanT5ft0UfL9ssv7Uwyp1BTI+ueg++27UqJWwvgP/+RpIZYITkZmCGNvbLzyildVB3dCRlCjhGQKFO6RgSk3JVAcg9n/Tn23nv3n2PBffXhh7Ldf3/oJ/VsHHOMbIuLZXlbJ/HOOyJ4iYlw3nm+z7EEZNs2Wd8kljjroWzySQFg7k3GCulqOOuO1Q2pzhMBqUpylvsKICfHvgmD8wVEa5g1S/ZPPdU+PnKk3dvLaW4sy3116qn+06NHjLCbQ8aaG2vgIMXmA8QKSf62kIaixnZGGGIJIyBRptFThd6QlRzlmfjGew1xpwvIypV2J9tTTrGPK2VbIU4KpOfl2RZF6+B5aywr5O23Y6c1i8VxD/ejingStObbW01hYVfCCEiUcRXVA6D6O1NA+vSx9ydOjN48AsGyPoYM2bPhoxUH+fprqax3Albqbv/+cOyxbZ97zjnS4qS62v6escLkw+NYPlgysupeK8Dd6LBAlCFojIBEmZRKEZCUYc4UEGvdc5fLbo/uVKwb6ymn7Nl917pB19ZKn6xo09wMzz0n+5dcIuLQFr16wYknyv7LL4d1amFh0p8H0gKkNzay4L6ids83xAYBCYhSKkkp9YBSqkApVaeUmq+UOjrAsQOVUm8qpcqVUpVKqfeVUrk+zvuDUuoDpdQOpZRWSt3Vwe8Sc7jdmuxGEZDMcc4UEGttbrfbdg85kR07ZC0N2N19ZdGnDxx4oOx//HHk5uWPTz6RoDjAFVcENuaii2T7+edQEGOeoOMuSmZphpiz2x/JR7tNYWFXIFAL5AXgBuAV4DrADXyqlDq4rUFKqR7AHGAqcC9wJ7A/8JVSKqvV6X8GpgBLAp18rLN9ZRNJiDnff4LzBETr3VfE+/HH6M2lPT76SObbsycccYTvc37xC9k6QUCeflq206fb/bra4+STJbGhpcUOvscKSsHgWyWTIbuylp+eNoWFXYF2BUQpNRk4F7hVa32r1voZ4CggH3igneFXASOBE7XWD2mtHwGmAwMRQfImV2vdD/CTzNj1yF9Qv2t/2BTnCcjmzbs/6TpZQD74QLYnnCApsb446STZrl8vDSKjxebN8Omnsv+rXwU+LinJrhX517+cV9PSHmfcms6yZEk1W3vXZrNWSBcgEAvkTKAJeNY6oLWuB54DDlNK9W9n7Hyt9S6rQmu9GvgSONv7RK31psCn3TUoXC4CUqXiyRrUjhM8CsybJ1vrhuxUAampsdNzfbmvLCZOtNOSP/oo/PPyx3PPibXUr59YFR3Bcndt3gz//W/o5xZO4uOhx5VDAMgpqmLd22atkFgnEAGZCKzWWrcuI10IKGCCr0FKKRcwHljk4+2FwGillMPDsuGlcq0ISHmy86wPkKpugFGjZPvjj85cXvWLL6C+Xm5QVidbX7hcthvLslgijXfw/LLL7MaJgTJ6NEybJvvPPBPSqUWE8+7PYE1CTwB+vDk/yrMxdJZABKQ/sN3HcevYAD/jsoGkNsYqz7WDxhOY9/sCHN14ozFfBMSpNSCWBXL44bItLZWutk7DagUybZrdONEfp58u22+/hcLCsE7LJx99JG5BpeDyy4O7xpVXynbWLLvNTKyQmqpImClWSL/8MtZ/UhXlGRk6QyACkgL4ypyv93rf3ziCHNstcHINSFWV3f31lFMkOA3w/ffRm5Mvampsa+Lcc9s//+ijIT1996r1SPKPf8j2uONg2LDgrnHGGZLWG4vBdIBzH+tFfnwaAAuv3hTdyRg6RSACUodYEq1J9nrf3ziCHBsQWuvMtl5ARWeuH26cXAMyb54EaV0uOOQQOOgg+7iT+Ogjqe1ISLCti7ZITrbdWO+9F965tWbFCjtWc+21wV+ndTA9lhosAqSmKVouHAbAgE0lrJtVGd0JGYImEAHZjm9Xk3XMX0Z6KWJ9+Bur8e3e6ha0NHvVgIxxnoBYLTYmThTr45BD5GenWSCvvy7b444LfKndM86Q7ezZUB7BOO5jj8l2zBiZb2ew3FibNsVeZTrAOU/kkBffA4DFVzt0pS9DuwQiIEuBsZ6aDm+meLbLfA3SWruBnwBfHZSmAOu01jHW1Sd0bP6hYVcNyKApzvPkzZkj2yOPlK0lIMuXi3vLCVRUSEEeBOa+sjjhBEhJgaYmePfd8MytNTt32hXk110nll1nGD3aTkv+6187d61okJqm0L+UAph+W8v4+T8mIysWCeSf8dtAArAr5KeUSgIuBeZqrQs8x4Yopcb6GHuQUmqi19gxSB1JjC3SGVo2fi3a2QKMmuYsAamqspewtQRkyhQJ/LrdsHBh9ObmzfvvywJSycltp++2pkcP+/zXXgvP3Frz9NPSgyszE2bODM01b75ZtvPmOc+1GAjnPZrNukQJri2/Ls/UhcQg7QqI1noBcrN/0NPO5Ergf8BQ4DavU18CVrUa/iSwEfhEKXWzUup64AvEdfWI94lKqYuUUn8EPP8tOFwp9UfPy9HZVMFQtFgEpCQ+mZSezmpJ9t134lePi4PDDpNjPXvCvvvKvpXeG20s99XJJ0tgvCOcf75s//c/aYMSThob7eD5FVdAWlpornv44XaH5IcfDs01I0lKiiLzRrFC+hdXMO/R0ijPyNBRAr1zzQQe82wfRyySE7XWbd5KtNZVwDTgO+AO4B7EJXaE1rp1L4PLPO//wfPzkZ6f7wHaSc6MPWrXiIBUZTuvFMZyX02aZGdfARx6qGy//jryc2rNzp1S/wEdc19ZHHecWANut50GHC7efFNEKi4OrrkmdNdVCm66Sfbfe89eVz2WmHFvFuvSMwHY+ocNuJtirLy+mxOQgGit67XWt2it+2utk7XWk7XWs1udM01rrXyM3aq1PktrnaG1Ttdan6K13qOawBrv57Up6G/oUFzbPAlog50nINaaGVbBmoXVEv2776K/JsVbb4mVlJ7edvGgP5KS4KyzZP/FF0M7N2/cbnjA0/DnjDOk1XwoOfNMuabW8Oijob12JHC5YOTDI3ADfetq+fQ3MdYlspvjLN9JNyKjQu7APcY6S0B27LBblrTOFDrySHnqbWyMvhvr3/+W7YwZEhAPBmsRpx9/hKVLQzOv1rz/vt2Q8ne/C/314+PhBk9XuX//OzrFkZ3lyCvS+WmIJ1nz+U1UFzRFd0KGgDECEgUqi1ro1SL1lX0PdJaAWP2V0tLs+IdFdrasNQ7RXRp22TJYvFj2L7ss+OscdBCM9aR9WO1FQonWcPfdsn/SSfbvLtRcdpn8berq4L77wvMZ4eboV3OpJo40dzMfnG7SemMFIyBRYO0c2/8z8khnZWBZXWKPPtp3V1tradhoCoh1sx8zxo7LBINStgC9+qrcgEPJhx+K2AHcb8OvnAAAGT9JREFUcUdor+1Nejr8/vey/89/SqPFWGOfwxLJP3IYAP0WFrDyo9at9wxOxAhIFNg6T+5UNcQxYB8/vcejQEuLbYH4iytYArJkCRQXR2Ze3tTVwSuvyP4vf7nnyoMdZeZMqWIvK4M33uj8/Cy0hnvukf3jj4fJk0N3bV9cfTUMGCDuRcvqiTUuemcg2+NTiAMWX7gWd7NJ63U6RkCiQPlyTwpvWiouVyfvgCFk4UJ7BUJ/AnLoobK0rdbRWZjptdfkZp+YCBdf3Pnr9ekDZ3sWFnj88dB1G/7sM7uW5v/+LzTXbIuUFPtzXngBVq8O/2eGmvQsF5l3jQZgcEUlsy7dGuUZGdrDCEgUaF4j5nnjgBAVBIQIqyp7n31g6FDf56SkyBM1SIA4kmgNf/+77J99NvTtG5rrWn2pliwJTXJASwv8wZOMfswxcHCb63aGjl/+EkaMkMyvSIhWODjhD1msGCYB9ZRX8ti6MMR+RUNIMQISBTKLpBdI2oQOVr+FEa3hnXdkf8aMts897TTZfv55ZNN5v//ezpYKZT3F5MlSaQ/wt791/novvihiBPDnP3f+eoGSkGC7r956y3l9ywLl5I9HUKySSMLNnBPX4G4xriynYgQkwhSuayLHk4E15JjW7cWix5IlkOdJfmlPQE46SYri6uoiuyqeVW19wAGhjynceKNs33sPfv45+OtUVdnWxwUX2MIUKc491872uuaa2OvUCzB0r3i42ePKKinnzQu2RXlGBn8YAYkwK94T66MF2Pc05wiIZX2MGiUurLbIyrKLDK1x4Wb1arv9+s03dz543poZM6RBIcD99wd/nb/8RWppUlI6d51gcbngiSdk/8cfY3PVQoAzH+zF2lzxUfZ6YwPL33FIB0/DbhgBiTDbvpL4R1FSKj37xEV5NoLWdjuPGTMCuzlbVdzvvhuZ7rwPPSTzHDGifQspGOLi7EK/116DtWs7fo3Nm20r6ZZbYPDg0M2vIxx8sL1eyG23wdYYjUXP+HoU2+NTSECz7vyfqdpuCgydhhGQCNO4Qu62NQOcE/+YO9fuo2Q1GWyPc86RLri1tfD22+GbG8iaF1Yr9FtukerrcHDhhZCbK0Ho22/v+PibbpKOuwMGwK23hn5+HeGhh6B3bxH3X//amWvZt0evwfEM+tfeNOCiV2M970xeg9sdg1+kC2MEJMKk7xALJHm8c9xX1rKokybZHXfbIzPTXgHwhRfCMq1d/OlPsnbHoEGhSd31R0KCHfR+5x2YPz/wse+9Z7vzHnggdB13gyUnx3Zlffyx3fol1ph6SQ8Kzh4FwLCtxbzyiy1RnpHBGyMgEWRnXhN9myQtcdA0ZwhITY3tvrJ6QwWK5Sb55htYty6k09rFqlXw0kuyf+edYvWEE+8g9I03ijXSHmVlUsgHMH26BM+dwFln2a7Ga6+FNWuiO59gueT1/qwb0Q+AIZ9t5JMbY7DhVxfFCEgEWfBPWXWtEcWkC3u2c3ZkeO01qK6WwrzzzuvY2KOPtutFHnmk7XOD5dZb5SY+alR4rQ8Ll8uOY3z/ffs9srSGX/0Ktm8Xq+Ppp0Mf4A8WpWQ+gweLq/Hss+WBIdZQCi5cPIq8dFkWKOGR1cz9e1mUZ2UAIyARZcen8o++IKMn6TnRD6Brbdc9nHNO4GuKW8TF2Z1gn39e1ugIJR9/DB99JPv33y8upkgwbRpcdJHs33Zb2x1uX3xRai5ARHTYsHDPrmNkZUmfr7g4WY74iitiMx6SkhHH9IX7sDUhlQQ0FdeuYMGrJjMr2hgBiSBpazzrPu/vjPWxPv9cXERgL0zUUS67TOIh9fX2qnuhoK4Orr9e9o8+WtbSiCR//avcfMvKpMLb10132TK46irZP/10uPzyPc9xAlOn2lbVa6/ZPbpijYFjEzjgy/GUuJJIpYWii5ax6JXKaE+rW2MEJEJsXtRA/0Yp2x5xljME5KGHZHv00bDffsFdo0cP+M1vZP+RR0Jnhfzxj5IZFh8vPaoi7Rbq0weefFL2P/nE3rcoKRHRqKuTzK1nn3WO68oX115rx6zuvDM87esjwdipyYz5YDzlrgTSdTM7L1rG3CeMOytaGAGJEIufk3/ktcQx+eLop/B+8YWsBw5SmNcZbr5ZrJDKytB0gv32WzumcscdsNdenb9mMJx7rh0Qv+EGmDdP9uvq4JRTpHI/JUUysDrq/os0Vjzk2GPl5yuusLsaxxr7nJTG8I8mstOVRAot1Px2Oe/fFGL/qSEgjIBEiNKPZQn47b0zSEqN7q/d7baL5g4/fM+VBztKdra93sVTT9kr8AVDYaHcuLWWtGJrnYto8eSTMG6cpBGfcYZYReefL2KilGSIBWu9RZrEREk1Pugg+f1efDH861/RnlVwjD8hlX3mTGRHQgqJaDL/9jMvHb6RlkazpnokMQISAYo3NzFkiyyekXlanyjPRm561rK1DzwQGtfL1VdLlXhLi9yYmoIoGm5qEvEoKJBFkl59NXKBc3/07AkffCAWVmEhTJhgdyF++GFZkzyWSE+XVvNTpsiDxJVXiksrkHRlpzHu8GSO+mkimzMkO2vIt/m83H85+csbozyz7oMRkAjw5R07SURTj4tj7s6J6lwKCuzMqXPOkafRUJCUJJlYSok43XVXx8a73RKs/uor+fn552XFQScwapQEn10uOw32xhvt32OskZEhK0pOny4/3323tIcpi8FQQp8xiZyzbT9W7TsIgGGl5SyZ8AMf/nYHOhbTzWIMIyARoO4DyQPdmtubjH5h6sMRAG63+L7Ly6FXLwlOh5KpU+1srvvug9dfD3xe11xj++T/9Kfw9LsKlg0bxOXn/ZT+1lv2crWxSI8esuTulVfKz++/D3vvbadNxxLJaS5+s3wkVTftRY2KI0M3kf7Eal7svYzln0RwvYHuiNa6y76A8oyMDB1NFr5SqecwR89hjv7i/pKozuX227UW77fWr70Wns9oaNB62jT5jKQkrT/7rO3z6+q0Pv98e15XXaW12x2euXUUt1t+TxkZMre4OK1vuknrHj3k55QUrf/9b+fMN1iefVbrtDT7bzBzptYFBdGeVXDkL6nXT/Vbsev/3Gzm6OfGrtIbvquN9tRiioyMDA2U6/buse2dEMsvJwjIc9lL9Rzm6FeTFuimhujdaf7xD/sGcc014f2skhKtR4+Wz0pI8C9WmzZpfcAB9ryuvVbrlpbwzi1QNm7U+rTT7Ln166f111/Le0uXap2ba7930klyfiyzcaPWRx5pf6eUFK1vvlnroqJoz6zjuN1af/yHnfqt+Hn2wxtf6adzV+oFL5Rrd6wrfgQIqYAAScADQAFQB8wHjg5w7EDgTaAcqATeB3L9nHsZsAqoB9YCVwfyGW18dlQF5PM/l+z6B/zpXcVRmUNLi9b33GPfGI49VuumpvB/7rZtWu+9t/25V1+tdXW1vNfUpPUTT9hP8kppff/9zniSX79e69/8Rqwna+4nn6z19u27n1daqvWMGfY5iYkigBs2RGfeoaClReunntI6J8f+XmlpWl9+udbff++Mv09HqKtq0S+duVW/7bKFZA5z9CvJC/WLJ+frjd/XRXuKjiXUAvIa0Ag8CFwJzPP8fHA743p4hGA7cAtwA5APbAayWp37K0B7xOYK4CXPzzcFMkc/nx81AVn7ba1+T32n5zBH/ztzyf+3d+5BcpVVAv+dnul5ZCaZvCZx8mJNkISsrAthiEZ2jeISRTekCrNVuqvilqQKyhV3XUULkNrSpVjWlfgIoiVKhS1cHu6iLKZggVhrEiQSCSRhCRp5xEwgmYQMk3lkpqfP/nHuzdy5093T3TPTPZk+v6pbt/t85+v7fafv/c793jowUPqn75VXVD/wgcHCYM0a1a6u0l2/vV31kksGrz9vnhVGS5cOyubOVX3kkdKlKRO9vaoPP6y6fr1qIjGYtvnzVe+5J3vBmU6r3nuv6oIFg3FEVNeutd/rOUPLp85O1ZtvVp0xYzBfoLp8ueoNN6ju2KGaSpU7lfnT3TGg93/ikN5dt3OII9nKVt1c/2v9wbt/r1u/cVw728+gTI0z+ToQUc09UkFELgKeAv5eVTcGsjpgL9Cmqn+eI+4XgVuAFar6TCBbFsS9WVW/EsjqgYPANlVdF4n/78BaYKGqduRMaObrn2hqamo6ceJEoVFHxZ7/7mLvFfto6eumU6pZsuUClq+ZUrLr799v8zHuuMP2pwBbcmTTJhstVUo6Omy00ubNw7dX/dCHbMe8efNKmyZVs9GTT8Ljj1tn8puRFTEWLbI0X3VVfsuyd3ebbb/9bTgYWW28ocEm7l16KbS22lL5pbb/aOjosP/tzjuHDxiYORNWrbKthVtbba7MwoU2Um2ikk4rO37YyYvfOMysF9pp0qFjzVMIbVMaSS1soH55A7NbG1j8/gb+6IIaqqom8DID48D06dPp6OjoUNXpufTycSC3Ap8DZqrqyYj8y8A/A/NV9XCWuDuBlKquiskfAc5S1WXB98uAh4E1qvpoRO9dWG3no6qa55ieIdcpmQM51Z3m2f/sYs9tr7PgN23UkqYfIbnxHay+Nud/UDSqcOSI7YR34AA8/bTNMN+zZ1CnpQU2brSVWMeL/n5ob7dlTA4etKXdX3zRhvPu2gWpVPa4tbVWAK1YYfNIliyB+fOtgJoxwwrhQuapqEJfny01cuTI4HHokKVr/37b8zw+ZFXERpFt2GC2Kmb+SSpls9I3bbLZ9PG5FcmkFbRhPs86yzZ9am62/Tuam22eRn29LX44UVC1//LHP4YtW+D55zPr1dbC4sVw9tm2vEtzsy0JE+axqclGfzU02Lm+vnzLv/T1pPnld9/k1XvaSe57g3m9XVmHpPaSoKOmlu6GWvpn1JJorqW+uZq65iR1zdXUz6mmsaWaafOqmTa7iikzEtQ3JUhOERKJM9PxjKUD+R9grqr+SUx+CfAYcJmqbskQLwF0A99X1c/Gwr4KXA80qmq3iFwPfA2YparHI3o1WJ/L11X1ugzXGMkzNDU1NVGMA/nbJUdY+vvXSKAIkECDgyGyKpRGUszmFDUM2vI1arklcS77EpntP4LZsxLGCxsWsiFiM49rakb+rULCVK1gjB8jMWWKOYWpU+1t/9gxW4AxH6qrLT+ZjjA9AwP5pyUkmbQJgnPmwNy5Y1s76Oszh3rkiDmrsCaYLyLmRKqq7K0+kRjMcxgePaKyYigk3sCAvTT095vTjNcsi71uNB/5pGs0zicetzHdz/J0B0v0JG/VLt7KSebTw2j8eBroI3H6SCGkTx+QRhiIfNfT3zmtF6LDzsMzn01n7q3ncOUXCttIJ18Hks+khBbgUAZ5WOvI1gAxE+t8z1Q7OQxI8NsHgvOpqPMAUNU+ETmW4xrjRuL1XlZyfGTFGIeoYwst/Bfz6U5X211UBlSt0Cq04BovurvtKIZcNZjR0N9vhfzRo1YrmUioWr7HK+8TheiLSrEvVWPBcZJsYzbbGJzom2SARfQwh17mcIpmTjGHXmbRRyMpppJiKv00ktmDJoA60tSVqxAI2PnKKDz8COTjQOqBTMVQbyQ8WzzyjFuPdcpnojfbNUasXlkNpSmXTjbqzp/GT55aePrNII2gQFokuB2EtJisO1HNiaoaDtU30pGsBRFminlQS0c0TZnSOXJY+Dn6ZlpdbTv01dYObXLJdb2RvsfbsDPpJ5ODR02NXT9smog3ORV6fVVbrDA8urutphIWpun04FvvwIClN5k0W0SPujqr9YQ7GA7tDp5Ye2Kk05bHnp6hb/bhOcxvNO3pdObPuWqOhTKWceK1xPBzPF/Z/qd8HU2msOL/6yq6aOQlGnkpl1paqU+nqE+nqB5IU0OaZDpNtQ5+rtEBqlBEg9YMHdqikYiGnQ4flhOAIXUPiZ0z6VzaOn7rAeXjQHqwmkScukh4tnjkGTfbNULdbNcYN77zy+nA+PRdOI4zmRAgGRyVRT5jJg5jTUxxQllblnjHsdpHtrjKYPPWYaBGRIYsih30gczKcQ3HcRynTOTjQHYDy0SkMSZfGZwzrgikqmlgD3BhhuCVwG9VNWwV3x2c47oXBmncjeM4jjOhyMeBPIDVzU5v2CkitcCngO2q2hbIFgVzPOJx3yki50fiLgXeB9wf0XsCq7FcE4t/NXASGDbKy3EcxykvI/aBqOpTInI/cKuIhKOmPgmcBVwZUd0MvIeh/Te3Y7PKfy4i/wakgH/Amqxui1yjR0RuBDaJyH3Ao8CfAX8DXKeqpZ0J6DiO44xIvmuLfwL4anCeATyHzf/YniuSqnaKyGrMWdyI1Xi2Ap9T1WMx3dtFpB/4PHA5NjP9WlUd40XHHcdxnLFgxImEZzLlWsrEcRznTGbMZqKfyYhIGpCmpqKmgjiO41QkHR0dAKqqOfvJJ7sDSWHNZm+OpJuF0PMUvJBjheL2Khy3WWG4vQqnGJtNA9KqmrObY1I7kNESrrU1UjXOMdxeheM2Kwy3V+GMp80m8OLLjuM4zkTGHYjjOI5TFO5AHMdxnKJwB+I4juMUhTsQx3EcpyjcgTiO4zhF4Q7EcRzHKQqfB+I4juMUhddAHMdxnKJwB+I4juMUhTsQx3EcpyjcgTiO4zhF4Q4kAyJSKyL/IiJtItIjIr8SkUvKna5yIyKrRUSzHMtiuqtEZJuIdIvIayLyTRGZUq60lwIRaRGRW0Rkq4h0BnZZnUV3rYj8RkR6ReRVEblJRIatfCoi00Xk+yJyVES6ROQJEfnTcc9MCcjXXiLycpZ77pYMupPZXq0isklEng/y9qqI/IeInJ1BN6/nb7RlXb47ElYadwFXABuB32Fb924Rkfeo6pNlTNdEYSOwKyZrCz8ED+zjwD5sC+MFwD8Ci4G/LFEay8FS4DrsnnkOWJVJSUQ+CDwIPAH8HXAe8BVgdvA91EsADwfhXweOAdcAvxCRFap6YNxyUhryslfALuy+i7I3+qUC7HUd8G7gfsxebwE+AzwjIhep6v9Bwc/fXYymrFNVPyIHcBGg2La7oawuMO7/ljt9ZbbN6sA260bQ+znwB6AxIvt0EPd95c7HONpnKjAr+LwuyO/qDHr7sAKxKiL7GjAAvC0i+6u4vYFm4A1gc7nzW0J7vQw8mMfvTXZ7rQJqYrK3Ab3AXRFZXs/fWJR13oQ1nI8A/cAPQoGq9gJ3AheLSEu5EjaREJGpWZpcpgF/gT2wJyNBm4GT2EM+KVHVTlU9lktHRJYDy4HvqepAJOh2rEn5iojsI1jN7qeRaxwF7gPWiUhyrNJeDvKxV5SguSVXM+hkt9cOVe2LyX6LvZCcCwU/f6Mu69yBDOd84IWY8QF2AgJMivbUUXI3tstjj4g8KiLnRcLOw5pGn45GCG783Zh9K5kw/3H7tGFvjefHdHdp8GoYYSf29j6s7XsScynQBXSJyAER2ZBBp+LsJSICzAXaA1Ehz9+oyzp3IMNpAQ5nkIeyeSVMy0SjD3gAuBa4HPgnrBq8TUTOCXTCt5ZsNqxk+0Fh9vF70XgOuAmrnV2FFZbfE5EvxfQq0V5/DczHallQ4vvLO9GHUw+cyiDvjYRXJKq6A9gREf1MRB7C3nZuwm7m0D7ZbFix9gsYyT5TYroVfy+q6trodxH5EbANuFFEvquq4V7fFWWvYOTjJswWdwfiQp6/UdvLayDD6QFqM8jrIuFOgKo+CzwGhEP/Qvtks2Gl268Q+/i9mIGg72gj5mzfFQmqGHuJyFuwEWdvAOtVNR0ElfT+cgcynMMMVgOjhLK2DGGVzkFgZvA5rP5ms2Gl268Q+/i9mJ2DwXlmRFYR9hKRJmAL0ASsUdXXIsElvb/cgQxnN7BMRBpj8pXB+dkSp+dMYDFwNPi8F0gBF0YVRKQG65TbXdqkTTjC/MftMw8br787prsi6CiNshIbUfO78UrkGcDi4Hw0Ipv09hKROuAh4Bzgw6q6P6ZSyPM36rLOHchwHgCS2LhpwIYPAp8CtgejZSoSEWnOILsYeC/wCEDQHv0Y8PHYjflxoBGbBFWxqOo+4AVgg4hURYKuBtLATyKyB7COzMtDgYjMBtYDP1XV/vFPcXkRkZnBBMGorA74AtAJRCe7TWp7BffLvViz3XpV/VVcp8Dnb9Rlne8HkgERuQ+b2HQbcAD4JNAKvFdVt5czbeVERJ4AurGO9Hbg7cAGoANoVdVXA70LAp292BjzBcDnga2qelkZkl4yROSG4OO5wMeAHwIvASdU9TuBzoeBn2Ez0e/F7PgZbG7INZHfqsI6SP8Ym1ndjs2sXgisUNXJ8Ead014iciVwPVbYvQzMwp7Hc4CrVfWOyG9NanuJyEZsBORDDI66Cjmpqg8Genk/f6Mu68o9u3IiHlgn0r9ibYS92Ljo95c7XeU+gM8CT2FLRPQDh7AHflEG3YuB7VhH3OvAt4CGcuehBDbSLMfLMb11wDPB/XUQGxJdneH3ZgSFQDs2D2IrcEG581kqewErggLzD9iIoTeBX2DNN5l+b9LaK8h3vvdXXs/faMs6r4E4juM4ReF9II7jOE5RuANxHMdxisIdiOM4jlMU7kAcx3GconAH4jiO4xSFOxDHcRynKNyBOI7jOEXhDsRxHMcpCncgjuM4TlG4A3Ecx3GK4v8BP0EFkGi7PTEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Start multiple series at 30, 60, 90, 120\n",
    "x = np.arange(200)\n",
    "stacked_y = np.zeros(len(x))\n",
    "for lambda_ in (30, 60, 90, 120):\n",
    "    stacked_y += plot_negative_binomial(x, lambda_, rho=2, color='b')\n",
    "plt.plot(x, stacked_y, \"-\", color='m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 417,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in real data\n",
    "def parse_hist(filename=\"reads.histo\"):\n",
    "    hist = [0] * 10002\n",
    "    with open(filename) as fp:\n",
    "        for row in fp:\n",
    "            cov, count = row.split()\n",
    "            cov, count = int(cov), int(count)\n",
    "            hist[cov] = count\n",
    "    return np.array(hist)\n",
    "\n",
    "hist = parse_hist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 298,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbd19f63850>]"
      ]
     },
     "execution_count": 298,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEXCAYAAAC9A7+nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZxjdZnv8c+TpJJautZeq/du6Ka7adZmEQUERXEUFUUdR8Rt0KuMIipedNSrOOpFcWYYBR1HRQd1HBAV5I64gIAiIns3NDS9QO9bdXdVakktqeR3/8hJkS6qupJUUklOvu/XK6+kT04lv5xOnjx5zm8x5xwiIlL5AqVugIiIFIYCuoiITyigi4j4hAK6iIhPKKCLiPiEArqIiE8ooIuI+ETJA7qZtZvZNWZ2j5n1mJkzs3Mm8XgBM/ugma01s14z22NmvzKzUwrYbBGRslPygA4cA1wFzAfWFeDxvgp823usjwP/BhwP3G9mxxbg8UVEypKVeqSomTUCYefcQTO7EPglcK5z7t48HisAdAN3OufemrF9NfAk8EXn3OcL03IRkfJS8gzdOdfjnDs40X5eKeVKM3vGzAa9Uso3zWxaxm4hoB7YN+rP93rX/QVqtohI2QmVugE5+D7wDuBG4DpgGfBhYJWZnedShszsQeA9ZvYX4I9AG/BFYA/wn6VpuohI8VVEQDezs4D3AG9xzv08Y/vDwH8D5wO/8Ta/C7gZ+HHGQ2wEznTO7ZmSBouIlEDJSy5ZegtwCLjPzGakL6Qy8ARwTsa+3cBTwDeBNwOXAbXAHWbWNqWtFhGZQhWRoZMqr7QBHePcPxPAzELA3cBdzrmPpe80s7uA9cAngM8Ut6kiIqVRKQE9QKoG/q5x7t/tXZ8NrAYuz7zTObfJzJ4BXla0FoqIlFilBPQtwLnAn5xzg0fYb7Z3HRzjvhoq5/WKiOSsUmrotwJh4FOj7zCziJk1ef/c6F2/fdQ+J5MawPR4MRspIlJKJR9YBGBmn/VuruSFronPA13Oueu9fb4LXArcQapOngSWA28DLnbO3eXt9zvgVcDPgbuAduAjpLLzk51zm6foZYmITKlyCejjNWKbc26xt48B/wv4AKnAP0gq6P8PcJ1z7oC3Xx1wJaksfYm335+AzzrnCjG1gIhIWSqLgC4iIpNXKTV0ERGZQEl6fZjZMKkvk+5SPL+ISIVqApLOuTFjd0lKLmaWBKy5uXnKn1tEpFJFo1EA55wbs7oyYYbuLTZxzzh3r3TObcijXd3Nzc3NXV1defypiEh1amlpIRqNjlvZyKXkch3w6Khtu8faUUREpl4uAf0+59xtRWuJiIhMSk69XMys0ZsAS0REykwuAf1HpHql9JvZ78zsuCK1SURE8pBNtj1Eai6VO4EDpBZcvpLUosunOuc2jv4DM5vobKe6t4iIFFhe3RbN7ATgEeAW59zFY9w/YUBvbm5GvVxERLLn9XKJOudaxro/r3q4c26tt2jEK8e5f8wnS/MCfsmz9Of6+9k7NMRL1R9eRHxgMkP/d5BaRagiJZzj7Mcf52WPP853d6v3pYhUvskE9KWMvyRc2Xuit5ddQ0MAfGDjRh5MjcASEalYEwZ0M5s5xrYzSa0g9NtiNGoq3Duqfn/FZk2TLiKVLZsa+s1mFgMeINXLZTWpOckPAF8oXtOK657OTgCmh0IcHB7msd5e4skkNQFNQCkilSmb6HUbMBP4BHADcBHwX8CpzrntRWxb0Qwnk/zJK7FcMX8+AHHn2NjfX8pmiYhMyoQB3Tn3Defc6c656c65GufcPOfc+yo1mEOqft6dSABwyZw5NAdTa0qv6+0tZbNERCalKusL93nZ+ZLaWhbV1rK6oQGAJ/v6StksEZFJqcqA/mwsBsCaxkYAjp82DVBAF5HKVpUBfefgIAALIhEAjktn6Cq5iEgFq+qAPn9UQN82OEh0eLhk7RIRmQwFdOA4r+QC8JTKLiJSoaouoPclEnR6WXg6oDeHQiPll/UK6CJSoaouoO/ysnN4IaADLK2tBWD7wMCUt0lEpBCqLqCnyy0GtIfDI9sXpgN6RsAXEakkVRvQ54TDhw3zX+hl68rQRaRSVW1Azyy3gDJ0Eal8CuiedIa+c3CQRB6rOImIlJoCuiedocedY583T7qISCVRQPcsyPi36ugiUokU0D2NoRCtodT08Kqji0glqqqAPphM0hGPAy8O6KCeLiJS2aoqoO/NqI1n9kFPU08XEalkVRXQD3nZOcCMmpoX3a8MXUQqWVUF9INeQA+Qmr9lNGXoIlLJqiqgH/Im5WoNhQiYvej+dIa+TRm6iFSg6groXobeNka5BV44UXpoeJgBb81REZFKUVUB/aCXoU8fJ6DPzej5skeDi0SkwlRVQB/J0Meon8PhPV92K6CLSIWpqoCePik6XoZeFwyODC7arROjIlJhqiqgp0+KjpehA8z1svRdCugiUmGqKqBPlKHDC3V0lVxEpNJUVUAfydCPFNC9DF0lFxGpNFUV0Ecy9COUXOYpQxeRClU1Ad05N2E/dMgouShDF5EKk1dAN7P/bWbOzJ4odIOKpTuRID1UKJuTosrQRaTS5BzQzWwO8Fmgr/DNKZ7MibmyOSnak0jQ49XcRUQqwfip6viuAR4h9WXQUtjmFM+hjOCczUlRSI0WbTxCNi8iUk5yytDN7DTgncDHi9Oc4kmfEA0CTcHguPvNyQjo6osuIpUk64BuZgZ8E/hP51zF1M7TMk+I2hgzLabVBALM8jJ41dFFpJLkUk94F7AKuHCiHc2sa4JdmnN43oKYaGKuTHMjEfbH4+rpIiIVJasM3cwaSdXOr3HO7Sluk4pjoom5Ms1TTxcRqUDZZuifBYaAf8lmZ+fcEU+Wehn8lGbp2Qz7T1NfdBGpRBMGdDNrB64APgfMzqg/1wJhM1sMRJ1znUVqY0FkMzFXmvqii0glyqbkMhsIA18Fns+4nA6s9G5fVawGFkpXevk5Zegi4lPZlFyeB940xvYvAQ3Ax4CNhWxUMXR7Af1IXRbTMjN059wRe8X40db+ftbHYpzS2MjsjG6cIlLeJgzozrkocNvo7WZ2BTDsnHvRfeUo6q0R2pRNycXL0AeSSbqGh7PK6ivdgaEhbty7l1s7Oni4p2dk+/mtrdy0ciWzFNhFyl7VDINMZ+jNOWTokBpc5PeAvndwkJc89hjbMkpMBjjgt52dnPboo9x5/PGsbGgoWRtFZGJ5z7bonDvHOXdiIRtTTN05ZOizwmHSYd/vJ0b7EgkuePJJtg0OUh8IcNncudxzwgn0nXUWP1m5kvpAgG2Dg7xm3Tr2+fxYiFS6qpg+1zlHNIcaesCM9io4MZpwjr97+mke7e0lANy8ahU3LF/OOa2t1AWDvGP2bO478UQaAgG2Dw7yxiefZCiZLHWzRWQcVRHQB5NJ4s4B0JzlZFvV0HXxY5s3c8fBgwBcv2wZF8yY8aJ9Tmlq4r9WrcKAv/b08K87d05xK0UkW1UR0NPlFsiu5AL+77p4x4EDfHPXLgA+uWABH5o3b9x93zBjBpd7939x61a2DwxMSRtFJDdVEdCjGVPnZlNyAX9n6B1DQ1z67LMAnNvSwjVLl074N1cvWcKccJhYMslVzz1X7CaKSB6qIqBnZuhZl1x8nKFfuWUL++NxGoNBfrBiBYEs+tk3h0Ijgf+W/fvZ0t9f7GaKSI6qI6BnZOjTqjxD/3M0yk379gFw7VFHsai2Nuu/vXjWLBZFIiSBr+/YUaQWiki+qiKgp0sujcFgVtkovJCh7xkaIumdUK10Cef48KZNAJzS2Mil7e05/X0oEOATCxYA8IM9e9SNUaTMVEVAT5dcsi23wAtT6A47R0fGeqSV7MY9e3iitxeAG5YtI5jHlAbva2+nLRRi0Dm+v6ciZ1IW8a3qCOg59EFPS2fo4I86enR4mM88/zwA7549m9OamvJ6nIZgkPfMmQPA9/bs8c2vFxE/qIqAnss8LmmtoRARL4P1Qx39y9u20RGP0xAI8JUserUcSbpU8/zAAH/oLOtZk0WqSlUE9FzmcUkzM9/0dNkci3GdNyDo04sWHfbrIx8rGxp4mZfhf1dlF5GyUR0BPY8MHfzT0+WTzz1H3DkWRiJ8fP78gjxmOku//cCBw3oRiUjpVEVAz2Uel0x+yNDv6+ritgMHAPjaUUdRl+MxGM+bZ84kYsagc9zuPb6IlFZVBPSRkkuVZejOOf73li0AvKSpibfNnFmwx24KhXjd9OkA3Lx/f8EeV0TyVx0BPd+Si5eh76rQDP2XBw7wkLdYxdeWLi34yktvnzULSM2ZfsgnXTtFKllVBPR8Sy7zKrjkMpxM8mlvzpULpk/nrJaWgj/H66ZPpyEQYNg5fqmyi0jJVUVAz2dgEbxQctkfjxOvsHnAb9y7l439/Rjwf5csKcpz1AeDI2UX1dFFSq86AvokT4o6qKhh7rFEgi9s3QrAu+fMYfW0aUV7rjd6c6j/vrOTvoxJ0ERk6vk+oB+2WlGeGTpU1onRa3fsYM/QEBEzrl68uKjP9TdtbQRJLah9lwYZiZSU7wN6fzJJOm/MteTSGAqNzM5YKXX0rf39XLN9OwBXzJ/PwhxmU8xHa00NL/fq8yq7iJSW7wN6dx6LW2SqtK6LV27ZwkAyydxwmM8uWjQlz/kGr+zy/w4eJKG5XURKxvcBPZrH8nOZKmlw0d2dnfzcy5KvPeoopuXxevPxBu/EaEc8zoPd3VPynCLyYr4P6JkZei5zuaSlp9HdWeYBPZ5Mcrk31/mZzc38nddHfCosqavjuIYGAH6lsotIyfg/oHsZupGa+jVX6RV9tpb5wsg37NrF07EYAeCbRx9d8EFEE3mDui+KlJzvA3rmoKJ8gtzSujoAnivjgL53cJDPe90UPzB3Lic2Nk55G9LdF5/t7+fZWGzKn19EqiCg5zuPS9pRXoa+c3CQwTIdXHTVc8/RnUgwPRTiy0UaRDSRNY2NtHvlKZVdRErD/wE9z3lc0tIZuqM8yy6Ziz5/ZelS2mpqStKOgBkXeGWX3xw6VJI2iFQ73wf0fOdxSZsXiVDjlWqe6+8vWLsKIeEc/7BxI5Ba9Pnvc1z0udDOb2sD4P5oVKNGRUpgwoBuZqeY2S/NbJuZ9ZvZXjP7jZm9dCoaOFn5zuOSFjRjsVd2Kbc6+rd27WJtXx+Q/6LPhfSKlhYCwJBz3NfVVdK2iFSjbDL0o4AQ8F3gw8C1wCzgj2b2qiK2rSDynccl01HpE6NllKFvisW4yptN8dL29rwXfS6k1poaTvfa8VuVXUSm3IRpq3PuZuDmzG1m9m3gOeCjwO+L07TCyHcel0xLvQx9yyQC+oPRKN/YtYvHe3tpDgY5ramJ98yZw8l59EiJJ5O8a8MG+pNJFkQifP2oo/JuV6Gd39bGX7q7FdBFSiCvGrpzLgZ0AIWfZLvAJltygcl3Xfz69u2c8fjj/HT/fjbEYvy1p4dv7trFmkcf5SObNuU8Ne8ntmwZGZH5gxUrJvXaCu3Vra1AqvvitjIrUYn4XdYB3cwazWyGmR1jZl8BVgN3F69phVGIkks6Q3+uvx+X41wlN+7Zwye90sgJDQ18belSPrdoEavq6wG4ftcuXr1uXdarIl2/cyff3LULgH9cuJBXegG0XJza2EiL9wXzO2XpIlMqlwz9B6Sy8g3AJ4B/B74y1o5m1nWkC9A86ZZnKTrJbovwQg29L5mkI4el1tb19vJBrxfKK1ta+OuaNXxy4UK+uGQJT556Kp9ZuBCAe7u6OP7hh/nJvn3jfmEMJBJ8aONGPrJ5MwCva2vjiyXqc34koUCA87wvGZVdRKZWLgH9auDVwPuAPwMRoDSdnnMwMrCoABk6wIYsR0HGk0nevWEDcedYVlfHL1evJhJ44XAHzPjS0qX84thjaQuFODQ8zDufeYazn3iCP3Z1jQT2pNdj5PTHHuPfd+8GUsPsf7pqVcl7tYwnXXa5u6uL4TIdjCXiR1mnrc65J4EnAczsx8AjwA+Bt4yx7xFr61OZpU92YBHAtFCIpbW1PDcwwNreXs7OYn3O63bu5IneXoxUnbtxnOd/08yZnNbUxEc3beLnBw5wfzTKy594grpAgJZQiL5EYuQ1BIHPL17MZxYtIlCmwRxe6I/eNTzMwz09nNE8ZT/IRKpavidF48DtwJvNrK6wTSoc51xBaugAJ3nLuD3e2zvhvv2JBF/fsQOAf5g3j5dNENDmRSLcuno1d59wAi/xuv31J5PsGRoaCeanNjbyp5NO4nOLF5d1MAdYWFvLCu8cgcouIlNnMt0j6khNYtgIlE8H7Qx9iQTpH/yT7QlyUmMjPz9wIKuA/sO9e9kfjxM249NenTwbr2ht5YGWFtb39fH8wADR4WHCgQBrGhtH6viV4vzWVjbEYvyus5MvlGGtX8SPJoxyZjbTOdcxalsT8FZgh3Nuf7EaN1ndk1zcItOJXoa+vq+PoWSScGDsHzfDySTXetn5u+bMGVkgI1tmxupp04q6sPNUeHVbG/+2axd/7e6mMx6ntURzzIhUk2yi3M1mNgA8AOwFFgDvBeYDby9i2yYtOsnFLTKlSy5x53i6r2/cKWp/1tHB8wMDGPDJBQsm9ZyV7OUtLYTNGHKOuzs7ecsULrghUq2yqaH/GKgHLge+DVwGrAXOdc7dUsS2TVohM/T2cJiZXpY5XtnFOTeyQPNFM2ey3KsjV6OGYJCzvHMHv+vsLHFrRKrDhAHdOXejc+4c59ws51yNc26mc+71zrn7pqKBk5E+IRoE6sYpkWTLzCY8MfqbQ4dY502WdVUVZ+dpr/Z6u/z20KGcB2SJSO58PX1uNGNxi0IsyZYO6H8ZZyHkdHZ+Xmsrp5TBZFmllu6+uH1wUKsYiUwBXwf0QvRBz5TOOB/p6XnRPCUPRKP8MRoF4FM59Gzxs+MbGpjtlalUdhEpPn8H9AL1QU87u7mZGV6A+kXHYR1/+KqXnZ/S2Mgrshh4VA3M7LCyi4gUl68DerQAMy1mCgUCvNlbDPlnGQH9r93d/OrgQSCVnReivOMX6bLLvV1dZbsmq4hf+DqgFzpDB3jLzJlAqo6+tb+fvkSCS555BoDjGhq40Av4kvIqb16XWDLJn72SlIgUh68DeiEWtxjt3JYWZnlll4vWr+dt69ezqb+fGjN+tHJl2U6YVSqzwmFO9k4mq+wiUly+DuiFWNxitFAgwPePOYYg8FhvL7/2gtSXlyzhhAof3VksqqOLTA1/B/QilFwALpgxg+96Qb09HOamFSu4Uv3Ox3W+V3ZZ29fH3iwX8hCR3JXP2mVFUIjFLcbz3vZ2/qatjdaamsPmOZcXe2lzMw2BAH3JJL/v7OSSOXNK3SQRX/J1JCrE4hZHMicSUTDPQjgQ4FytYiRSdL6ORoUeWCT5S5ddft/ZSVLTAIgUha8DerRINXTJXbo/+v54nLVZzCkvIrnzbUBPOkdPEXq5SH6Orqtjibc2q8ouIsXh24DeW8Cpc2XyzGxk8WjN6yJSHL4N6JmLW6jkUh7SZZf7o1F6M/5/RKQwfBvQMxe3UMmlPLyitZUgqVWf7u3qKnVzRHzHvwFdGXrZaQ6FeIk3T7zKLiKF59uAni651Jipr3gZOV/TAIgUjW8jXeY8LprOtnykA/rG/n629veXuDUi/uLfgK4+6GVpTWMjbd45DZVdRArLtwG9mPO4SP6CZpynaQBEisK3Ab3Y87hI/tJll7s7OxnWKkYiBePfgK4MvWylBxhFEwke6ukpcWtE/MO3AV3zuJSv+bW1rKqvB+DX3lqsIjJ5vg3oIyUXZehl6fXTpwNw24EDJW6JiH/4N6Cr5FLW3uQttr0+FmNTLFbi1oj4g28Dukou5e3Uxkbaw2EAfqksXaQgfBvQi7FAtBROwIwLZ8wAVHYRKZQJA7qZnWpmN5jZ02bWZ2bbzey/zezoqWhgvjSwqPy9yQvof+nuZo8WjxaZtGwy9KuANwN3AR8F/gM4B3jczFYWr2mTM1JyUYZets5paRkZJ3C7snSRScsmoP8LsMg5d7lz7nvOuS8BZwE1pIJ92Uk4R583YEUll/JVEwhwgdfbRXV0kcmbMKA75x5wzg2N2rYJWA+UZYbeo6lzK0a6t8sfurroisdL3BqRypbXSVFLTV84GyjLtCqq5ecqxmva2qgNBBh2jl9rbheRScm3l8vFwDzglrHuNLOuI12A5nwbnI3MxS00l0t5awgGeZU3FYDKLiKTk3NAN7MVwA3A/cCPCt6iAshcfq5RGXrZS/d2ufPgQfoz/u9EJDc5BXQzmwP8D9AJvNU5N+ZUec65liNdgOjkmz6+dA+XiFYrqgivnz6dANCXTHKX5kgXyVvW0c7MmoE7SZVLznfO7S1aqyapW10WK8qMcJizmlNVOA0yEslfVgHdzGqBO4DlwAXOuWeL2qpJ6vICeqsCesVI93b51cGDmiNdJE/ZjBQNAjcDZ5AqszxY9FZNUqcX0FsU0CtGehqAA/E493Z1lbg1IpUpmwz9n4E3kCq3tJnZOzMuFxa3eflRhl55FtXW8tKmJgB+tG9fiVsjUpmyiXgnetev9y6ZtgG3FbRFBaAMvTJdMns2D3R38/OODr61fDkN6nIqkpNsRoqe45yzcS6Lp6CNORvJ0GtqStwSycXbZs0ibEZfMqmToyJ58GWfPmXolamtpmZkbpcf7S3bTlQiZcuXAb1LAb1iXTJ7NgC/7+zUlLoiOfJ1QNdJ0crz2unTaQuFSAI/3b+/1M0RqSi+DOid3qx9ytArTzgQ4G9nzQLgJpVdRHLiu4DunFOGXuHSZZe1fX2s6+0tcWtEKofvAnpvIkF6eidl6JXpJU1NLK+rA+A7u3eXuDUilcN3Ab0rY+pcZeiVycz44Ny5QGqQUeaCJSIyPl8HdGXoles9c+ZQFwjQk0jwE40cFcmK7wJ6Z+biFgroFau1poa3eydHv7V7N865ErdIpPz5LqCnM/RpwSA1mgu9ol3mlV2e7Ovjge7uErdGpPz5LuJplKh/nNLUxKmNjQB8a9euErdGpPz5LqCry6K/pLP0Wzs62D80VOLWiJQ33wZ0Zej+8LezZtEaCjHkHN9WF0aRI/JdQE+PElWG7g91weBIF8Zv7NxJr7owiozLdwFdGbr/XDF/PrWBAIeGh/nunj2lbo5I2fJdQNdJUf+ZFQ5zaXs7AP+8YweDWnNUZEy+C+g6KepPVy5YQMiMXUNDmitdZBy+DejK0P1lUW0t7/AGGn1txw4SGmgk8iK+C+gqufjXpxYuxIBN/f38THOli7yI7wL6Qa+Xy3StJ+o7KxsaeNOMGQB8butW4qqlixzGVwE9lkgQ8z7kMxXQfelLS5YQADb392tqXZFRfBXQO7zsHGBmOFzClkixrGxoGOnxcvW2bXRl/J+LVDt/BfSMoeHK0P3r6sWLmRYMciAe53Nbt5a6OSJlw1cB/YCXrdWY0RgMlrg1UixzIhG+sHgxkJq067GentI2SKRM+Cqgp0suM2tqMLMSt0aK6fJ58zi2vp4k8A+bNpFUN0YR/wZ08beaQIAbli8H4MHubn6gwUYi/gro6ZLLDAX0qvDylhbeOXs2AFdt2XLYORSRauSrgD6SoauHS9W4dulSmoNBDg4P8/5nn9VSdVLVsgroZtZuZteY2T1m1mNmzszOKXLbcpbO0FRyqR5zIhGuX7YMgNsPHuR7mo1Rqli2GfoxwFXAfGBd8ZozOR0quVSli2fP5m9nzgTgis2b2RiLlbhFIqWRbUB/FJjhnFsGXFvE9kzKAZ0UrUpmxreXL2d+JEIsmeSdzzyjaQGkKmUV0J1zPc65g8VuzGSpl0v1aq2p4aYVKzDg4Z4ePr5lS6mbJDLlfHNSNJ5Mjsy0qJOi1enc1lb+ceFCAK7ftYvvaa4XqTJFmWPWzLom2KW50M95KGOtSdXQq9cXlyxhXV8fdxw8yGWbNrGivp4zW1pK3SyRKeGbDF3zuAhAwIwfr1zJqvp64s5x0fr1bB8YKHWzRKZEUQK6c67lSBcgWujnTNfPDWjT4hZVrSkU4vbVq2kNhdgfj/OadetG5skX8TP/ZOjeB7Y1FCIU8M3LkjwdXV/PL449lrAZz8RivHbdOk21K77nm8inLosy2jmtrfzXqlUY8FBPD+etXatMXXzNNwF9j1dDn60eLpLhopkz+fHKlQSBR3t7OeOxx9ikgUfiU1kXm83ss97Nld71JWZ2JtDlnLu+4C3L0Q7vxNeCSKTELZFy847Zs4kEAlz89NNs6u/nJY89xi9Xr+Zs9X4Rn8nl7OE/jfr3+7zrbUDJA/r2wUEAFtbWlrglUo4umjmTuSeeyBufeoqOeJzz1q7lW8uW8fft7Zo7X3wj65KLc87GuSwuYvuylu6atlAZuozjjOZmHjz5ZFZ6XRrfv3EjFz/zDN0ZYxhEKpkvauhJ59jpZegLlKHLESytq+OBk07izTNmAPDT/ftZ8+ijPNTdXeKWiUyeLwJ6RzzOoDcPtjJ0mUhLTQ23HnssNyxbRsSMzV5d/QPPPssBLZIhFcwXAT1zJKBOiko2zIzL5s3jwZNP5riGBhzw3T17WP7QQ1y/c6dma5SK5IuAvsMrtzQGgzRrlKjk4MTGRh5bs4ZvHH00zcEgncPDfGTzZlY9/DDf2b2bvkSi1E0UyZovAnrmCVH1WJBchQIBPjJ/PhtPP533t7cTADb39/PBjRuZ/5e/cOXmzTzV26vl7aTs+SKg79AJUSmAWeEw/3HMMTx16qm8d84cImZ0DQ/zzzt3ctwjj3Dsww/zheef5+m+vlI3VWRMvgjoI33QVT+XAljZ0MCNK1aw44wz+PKSJSzxEoVnYjGu3raNYx9+mNUPPcTVW7fyRE8PSWXuUiZ8UXDerlGiUgQzw2H+cdEiPr1wIY/29HBLRwe37N/PtsFB1sdirN+6lS9s3cqMmhpe2dLCea2tnNfayuK6ulI3XaqULwL6Do0SlSIyM05pauKUpia+unQpD/f08LOODm7t6GDrwAAH4nFu7ujg5o4OABZFIpzV0sJZzc2c0tjIqvp6aoPBEr8KqQZWihM9ZtbV3Nzc3NU10cJGE+tPJGj4059wwB9OOIFzW1sn30CRLM4kwoYAAAzJSURBVDjn2NLfz12dndzV2ckfurpGlkHMFARW1NdzbEMDxzY0cHZzM2saG2lUjyzJUUtLC9FoNOqtK/EiFf+Oeqqvj/RX0qqGhpK2RaqLmXF0fT1H19fzwXnzSDjHE729/Kmriz9Go/w5GmV/PE4CUiWaWAy8LB5gdk0Ny+rrWVZXN3JZ2dDA8ro6ajSnv+Sh4gP6E729ALSHw5o6V0oqaMaaxkbWNDZyxYIFOOfYOzTE2t5e1vX1sSEW45GeHp70esnsi8fZF41yf/TwBbzCZqxqaGB1QwNH19VxVG0tR3vBXuMs5Egq/t3xuBfQT5o2rcQtETmcmdEeidAeifCa6dNHth+Kx3k2FmNTf3/qEoux0bvdm0gw5GX66WQl0/xIhGPr61lRX88xGdft4bDGYIgCushUa6up4YzmZs5obj5su3OOHYODrO3tZW1vLxtiMbYMDLClv39kicWdg4PsHBzkt52dh/1tYzD4QoCvqxsJ9Mvq6nRCtopUdEBPOMc6L6CfqIAuFc7MWFhby8LaWl7vzQaZ1hmP83Qsxvq+Pp7u6+PZ/n42xGJsGxjAAT2JBI/09PBIT8/hjwksrq3lmPp6ltfVsby+niW1tcwNhzmqrk4nZn2mov83N8VixLxJlE5qbCxxa0SKp7Wmhpc1N/OyUVl9fyLBpv5+no3F2BCLjVxviMXoSyZxwPMDAzw/MMBvxnjc9nB4JJs/pq4uFfjr61lcW0tQJZyKU9EBPV1jbAoGR0bziVSTumCQ46dN4/hRv1Cdc+weGhoJ8hu9Ov3GWIwdg4MMed2V9wwNsWdoiHtGdSEOm3FUXd1IVr+4tpYFkUjqF0QkQksopJp9GarogP6Q9/PyhGnTCOjNJTLCzJgXiTAvEuGVo8ZmOOc4GI+z0SvbZGb3WwYGGHaOIed4JhbjmVgMDh580eNPCwZTAd4L8pnBfkEkwvxIRLX7EqjYgD6UTPKTffsAXvSGFZHxmRkzwmFmhMO8dFQJZziZZOvAAM962fxGrxfO9sFBtg8MjCwk05tIvBDwx9EaCtEeDjM3EqE9HB7pWjwn4zI7HKZN2X7BVGxA/0VHB/vjcYLApe3tpW6OiC+EAoGRwVKvy+hqCanM/kA8PhLctw8OssO73j4wwI7BQfYMDY0M9OscHqZzeJinjxD0AWrMmB0OM7um5kXBfvS/G4NBBf8jqLiAvq63l+f6+7l2xw4A3jhjBvM0KZdI0ZkZM8NhZobDrBmnE8JQMsmuwUF2DQ6y26vP7/Zu7xsaYq933RGPjwT+uLcmcHpd4COpCwQOC/aza2qYXlNDayhEa/o6FKKtpoZ54TDTa2qq6gug4gL65Zs2cV/GyLoPzZ1bwtaISKZwIMCSujqWTDDj5HAySUc8PhLk9w4NsS8eH7m9N+MLIHN+nP5kcqTXTjbqAgFaQiGagkGaQiGaM25nBv+xbjeHQoQrbAqGigvocyMRpgWDzAuHeVVbG69Q/Vyk4oQCgZFRtBMZTCbZP0agT38JdMbjHPLKO53xONGMZQP7k0n6h4bYk2c7018IzaEQLd6lORg8bNvo68x9Gqa4RFRxsy3Gk0lNXCQi40o4R9fwMNsHBtg5OEh3IkH38DDR4eGR213eJfOL4NDwMPECx8Mg0OwF+wCpXzBzwmHmRSLctGJFzsHed7MtKpiLyJEEzZju1dZzGXDonCOWTNIZj9PpfQGkA380kUhdp/89+tq7f8Ab6JiWAA55Xxxpz8RizCnS3DsVF9BFRIrBzGgIBmkIBpmf52MMJpNjBv1oIkHSOQaSSfYODVFTpDKMArqISIFEAgFmhcPMKtFU3qpfiIj4RFYB3cwiZvZVM9ttZv1m9qCZvbLYjRMRkexlm6H/EPgY8GPgo0ASuNPMzihSu0REJEcT1tDN7DTg7cDHnHPXedtuAp4CvgqcXdQWiohIVrLJ0N8CxIHvpTc45waA7wNnmpkmUhERKQPZ9HI5CdjgnBu9wOFDpBZEOREOH4hlZhONGGqe4H4REclRNgG9Hdg1xvZ0EM9rMpVoNEpLy5iDnUREZAzR1DxWTePdn01ArwPGmgZtIOP+w4w3LDXNzIaBQDQa7R5nl3QGHx3nfhmbjlt+dNzyo+OWn8kctyZSnVLGlE1A7wfGmkGnNuP+nDjnjvi86ZLNRF8Mcjgdt/zouOVHxy0/xTxu2ZwU3UOq7DJaetvuwjVHRETylU1AfwJYYWbTRm0/3bteW9gmiYhIPrIJ6LcCNcCl6Q1mFgHeC/zZOacMXUSkDExYQ3fO/dXMfgZ8zetzvgV4N7AIeE9xmyciItnKdrbFdwH/5F23AuuA1zrn/lyshomISG5KsmLRRHT2PD86bvnRccuPjlt+St3LRUREKkBZZugiIpI7ZegiIj6hgC4i4hMK6CIiPqGALiLiE2UT0LVu6ZGZ2Tlm5sa5rBi170vN7H4zi5nZXjP7NzOrL1Xbp5KZtZvZNWZ2j5n1eMfnnHH2fYOZPWZmA2a23cw+b2YvGpthZi1m9h9m1mFmfWb2BzM7segvZgple9zMbOs478FrxtjX18fNzE41sxvM7Gnv9W03s/82s6PH2Derz+Rk42C2A4umwg+Bi4DrgM2kRqHeaWYvd879pYTtKjfXAY+O2jYy/YL3gbkbWA98HJgPXAksBV4/RW0spWOAq0i9h9YBLx1rJzP7G+A24A/AR4DjgP8DzPD+nd4vAPyPd//XgYPAZcC9ZrbGObelaK9kamV13DyPknofZnoq8x9VctyuAl4G/IzUMZsDfBh43MxOc849Azl/Jn/IZOKgc67kF+A0wAFXZGyr9V7QH0vdvnK4AOd4x+jCCfb7NbATmJax7VLvb19R6tcxBcepEZju3b7Qe93njLHfelKBKZix7UtAAliWse1to487MBPoBG4q9estwXHbCtyWxeP5/riR+tILj9q2jNRaET/M2JbVZ7IQcbBcSi5atzQHZtY4TmmgCXgVqQ9M5pKBNwG9pD5kvuac63HOHTzSPma2ClgFfMc5l8i461ukypAXZWx7C6lfQLdnPEcHcAtwoZnVFKrtpZTNccvklQaOVMbz/XFzzj3gnBsatW0TqWRhJeT8mZx0HCyXgJ7NuqWS8iOgG+g3s9+Z2XEZ9x1Hqoz2SOYfeG+6J0gdZ3nhOIw+TrtJZVInjdr3UeelSxkeIpXVvqheWgVeDfQBfWa2xcw+MMY+VXnczMyA2cABb1Mun8lJx8FyCejtjFpo2jOpdUt9ZojUVMYfBd4IXE3qJ9r9Zrbc2yf9DT7esdRxTMnlOOm9ebh1wOdJ/Yp5P6nA9R0z+9So/ar1uF0MzCP1SwSm+L1WLidFc163tNo45x4AHsjY9Cszu4PUN//nSb2R0sdpvGNZ9cfRM9Fxqh+1r96bHufcGzL/bWY/AO4HPmdm33bOpdfJrLrj5vU2u4HU8fiRtzmXz+Skj1m5ZOgFX7e0Gjjn1gJ3AeluTenjNN6x1HFMyeU46b15BN45iOtIfQmekXFXVR03M5tDqldPJ/BW51x6Iecpfa+VS0DXuqX52wG0ebfTP83GO5Y6jim5HCe9Nye2w7tuy9hWNcfNzJqBO4Fm4Hzn3N6Mu6f0vVYuAV3rluZvKdDh3X4KGAZOydzBzMKkTqg8MbVNK1vp4zD6OM0l1Uf4iVH7rvFOdmU6nVQvhc3FamQFWepdd2Rsq4rjZma1wB3AcuAC59yzo3bJ5TM56ThYLgFd65ZOwMxmjrHtTOBc4LcAXv3yLuCSUW+KS4BppAZAVD3n3HpgA/ABMwtm3PUhIAn8PGPbraRORr0xvcHMZgBvBW53zsWL3+LyYGZt3oChzG21wCeBHiBz4Ivvj5v33rmZVKnprc65B0fvk+NnctJxsGzmQzezW0gNaPhXXli39FTgXKel7jCzPwAxUidGDwCrgQ8AUeBU59x2b7+TvX2eItWfdT7wCeAe59xrS9D0KWdmn/VurgTeAdwIPA90Oeeu9/a5APgVqZGiN5M6nh8m1Tf9sozHCpI6yXUsqRGPB0iNeFwArHHO+SLThImPm5m9B/gMqcCzFZhO6nO6HPiQc+7fMx7L98fNzK4j1evsDl7o1ZLW65y7zdsv68/kpONgqUdbjRoRdS2pOtIAqb6X55W6XeVyAS4H/kpqCHUc2EXqA7dwjH3PBP5M6iTKPuAbQEOpX8MUHis3zmXrqP0uBB733m87SHUFDY3xeK3eB/EAqf7X9wAnl/p1TvVxA9Z4wWsnqd4Y3cC9pEoNYz2er4+b99qzfa9l9ZmcbBwsmwxdREQmp1xq6CIiMkkK6CIiPqGALiLiEwroIiI+oYAuIuITCugiIj6hgC4i4hMK6CIiPqGALiLiEwroIiI+8f8B4wEfU95zWYYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(5, 201, dtype=np.int)\n",
    "y = hist[x] * x\n",
    "plt.plot(x, y, \"c-\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "46605857762\n",
      "29574390523\n"
     ]
    }
   ],
   "source": [
    "# Decompose the truncated histogram\n",
    "X = np.arange(1, len(hist), dtype=np.int)\n",
    "Y = hist[X] * X\n",
    "total = Y.sum()\n",
    "print(total)\n",
    "print(y.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 340,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import entropy\n",
    "\n",
    "coverage_range = np.arange(5, 201, dtype=np.int)\n",
    "P = Y[coverage_range]\n",
    "\n",
    "def generative_model(G, lambda_, rho):\n",
    "    stacked = np.zeros(len(coverage_range), dtype=np.float64)\n",
    "    n = lambda_ / (rho - 1)\n",
    "    p = 1 / rho\n",
    "    for i, g in enumerate(G):\n",
    "        if 1 < g < 9e9:\n",
    "            stacked += g * nbinom.pmf(coverage_range, n * (i + 1), p)\n",
    "    stacked *= coverage_range\n",
    "    return stacked\n",
    "\n",
    "def func(lambda_, rho, G):\n",
    "    stacked = generative_model(G, lambda_, rho)\n",
    "    return np.sum((P - stacked) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 321,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbd1a210fd0>]"
      ]
     },
     "execution_count": 321,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEXCAYAAAC9A7+nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZhcZZX48e+pvXpf0+nu7BtJIBIIEEAUEBVxEBkQnZ8LoqM4ODiiqKCDCiMq4oYLKrghOAoCgiKiI4ggOwkEyL6nO+l936prfX9/3FvVlU53eq2upc/neeqpzr23qk7dVJ96+9x3EWMMSimlsp8j3QEopZSaHprQlVIqR2hCV0qpHKEJXSmlcoQmdKWUyhGa0JVSKkdoQldKqRyR9oQuItUicpOIPC4ivSJiROSsKTyfQ0T+Q0ReEZE+EWkUkT+KyEnTGLZSSmWctCd04BjgGmAe8Oo0PN83gB/bz/Vp4HvA64CnROTYaXh+pZTKSJLukaIiUgh4jDHtInIh8ABwtjHmH5N4LgfQAzxijLkkaftxwGvA/xhjvjw9kSulVGZJewvdGNNrjGkf6zi7lPIZEdkmIkG7lPIDESlIOswF5AHNwx7eZN8HpilspZTKOK50BzABPwfeC/wCuAVYDlwJrBaRNxtLSESeAy4TkWeBJ4Ey4H+ARuBX6QldKaVSLysSuoi8AbgMeJcx5v6k7S8CdwPnAn+xN18K3AP8OukpdgJnGGMaZyRgpZRKg7SXXMbpXUAH8ISIVMRvWC3wKHBW0rE9wGbgB8BFwMcBH/CQiJTNaNRKKTWDsqKFjlVeKQNaR9lfCSAiLuAx4FFjzKfiO0XkUWALcDXw36kNVSml0iNbEroDqwZ+6Sj7G+z7NwLHAf+VvNMYs0tEtgGvT1mESimVZtmS0PcAZwP/NMYEj3JclX3vHGGfm+x5v0opNWHZUkO/D/AA1w7fISJeESmy/7nTvv+3YceciDWA6eVUBqmUUumU9oFFACJynf3jKoa6Ju4DuowxP7SP+SnwEeAhrDp5DFgBvBt4nzHmUfu4/wPeAtwPPApUA5/Aap2faIzZPUNvSymlZlSmJPTRgjhgjFlkHyPAx4DLsRJ/ECvpPwzcYoxps4/zA5/BaqUvto/7J3CdMWY6phZQSqmMlBEJXSml1NRlSw1dKaXUGNLS60NEIlhfJj3peH2llMpSRUDMGDNi7k5LyUVEYoAUFxfP+GsrpVS26u7uBjDGmBGrK2O20O3FJh4fZfcqY8z2ScTVU1xcXNzV1TWJhyql1OxUUlJCd3f3qJWNiZRcbgE2DtvWMNKBSimlZt5EEvoTxpgHUxaJUkqpKZlQLxcRKbQnwFJKKZVhJpLQ78LqlRIQkf8TkTUpikkppdQkjKe1HcKaS+URoA1rweXPYC26fLIxZufwB4jIWFc7tXuLUkpNs0l1WxSR44ENwO+MMe8bYf+YCb24uBjt5aKUUuNn93LpNsaUjLR/UvVwY8wr9qIR54yyf8QXi7MTftpa6aFYjL91dvLG4mIKXXpJQCmVG6Yy9L8eaxWhrPPdgwc5/7XXeOOmTfRHo+kORymlpsVUEvoSRl8SLqPd2dQEwKa+Pt6/dQtNTXcTCOxLc1RKKTU1YyZ0EakcYdsZWCsI/TUVQaVSTyTCjoGBxL8fbO/g7u03sXnzO9GZJ5VS2Ww8BeR7RGQAeAarl8txWHOStwHXpy601Hiiq4so4BKhXPppjuWxg2NY2/87+vpeorBwXbpDVEqpSRlPyeVBoBK4GrgVuBj4DXCyMaYuhbGlxGOdnQCcXJDHsbHnANjNCgCamu5KW1xKKTVVYyZ0Y8z3jTHrjTHlxhi3MabWGPPhbEzmAI/ZXSXPyA+wjF0A1LlPAaCl5bfEYpG0xaaUUlMxqxa46AiH2dzfD8Apzj0sZQ8Ae8KFhHATDrfQ1fVYOkNUSqlJm1UJfd/gYOLn+ZENLMNaLzoKtPjPA6C7+9l0hKaUUlM2qxJ6vZ3Q8xwOHH3PU0Yn5Y4QAHXu9QD0929OW3xKKTUVsyuhB4MALPB6GRiwEveaPAFgtywHoL//tfQEp5RSUzSrEnqdndBr3BGMsX4+obACgJ0R6z4Q2E00GkhPgEopNQWzKqHHSy5VYnVddLlKeV3RXAD2hd32UTEGBralIzyllJqS2ZXQ7RZ6hWkEoKDgeBb5fAA0hCIYVzWgdXSlVHaalQl9jrGWQvX5FiUSOkCv/zRA6+hKqew0axJ6JBbjkJ3Qy6P7AfB4aqn1ehMnocNzAgB9fZrQlVLZZ9Yk9MZQiJj9c2nEGiHq9dbidjiY5/UC0O7Sni5Kqew1axJ6vNwCUBq2auReby0AC+2ySzPWBdJQqIFodACllMomsyah19k9XMpcTrymGzgyoR+KFiSOHxzcP7MBKqXUFM2ahB5vode6h7Z5PFZCj18YPRh2IGKVXwKBvTMboFJKTdGsS+g1Lmuov4gLj2cOAAvtGvqBwUF8vkUADA7qCkZKqewyaxJ6vIdLlcOabdHjqUbEevvxksvBYBC3bykAg4PaQldKZZdZk9DbwmEAirHmQ4/Xz2Go5BIFut2rAXSNUaVU1pk1Cb3dTuiFpg0Aj6cmsW++XXIBaHVqC10plZ1mT0KPWCsRFcSagMNb6D6nk7keDwAtjnmAVUPXRaOVUtlkViR0Y0yihZ4fOQgcntBhqJXeasoBiEb7CIfbZjBKpZSamlmR0HujUcJ2azsvvB8Y6rIYFx8t2mIKE9u0p4tSKpvMioQeb50D5NnzuAxvodfaCb0hZHC5SgHti66Uyi6zLqEXEx8lWnPYMbV2Df1QMIjPtwTQC6NKqewyKxJ6vMuiC8jDmqPF7a487Jh4yeVQKJQ0uOjAjMWolFJTNSsSeryHS5nLYK0g6sDlKjnsmHjJpTkUwuldCEAwWDeDUSql1NTMjoRut9BLnVEA3O6yxCjRuHhCN0CPK94XXVvoSqnsMSsSerzkUuKwhv+7XOVHHFObNLioPdEXvU77oiulssasSOjxFnqJBABwuyuOOCbf6aTE5QKgTaoAiMX6iUQ6ZihKpZSamkkldBH5nIgYEdk03QGlQjyhF9ELgNt9ZAsdhnq6tJqh+rqWXZRS2WLCCV1E5gLXAf3TH05qxC+KxrssjtRCh6GeLo1hJw5HHmCVXZRSKhtMpoV+E7DBvmWFtsTEXO3AUVroh3VdjPd00Ra6Uio7TCihi8gpwPuBT6cmnNSIl1wKYtbcLKO10BMJPRjE610AaAtdKZU9XOM9UEQE+AHwK2PMJuufox7bNcbTFY/3dadDYurcaCMwegt9XlJC95VaLXStoSulssW4EzpwKbAauDBFsaREIBplIBYDID92CJhYC10HFymlssW4ErqIFGLVzm8yxjSOdbwxpuRo++0W/Iy00pPncSmINQMj90OHoV4uQWPody0GtIWulMoe462hXweEgO+kMJaUiPdwgfH3cgHosAcXhcMtRKOBFEaolFLTY8yELiLVwFXArUCViCwSkUWAD/DY/y5NaZRT0JncQqcPGL2GXu5247WvDbTJ0ORdwWB9CiNUSqnpMZ4WehXgAb4B7Eu6rQdW2T9fk6oAp6o7as3fUuAwOIkBkpjvfDgRoSY+SVeskPjpmS09XYKxGE90dbGpt5eYTnmgVNYZTw19H/CvI2y/EcgHPgXsnM6gplO3XXIpdMQgBi5XCQ7H6G97ntfLvsFBGkJhjvfWEgzW53RfdGMMf2hr457WVh5ub6fX/gKscrv5/vLlvHvOnDRHqJQarzETujGmG3hw+HYRuQqIGGOO2JdJhhK6dT9a/TwuuaeLz7eQYLA+Zy+MGmP41O7dfO/QoSP2NYfDvGfrVjb19fHVxYs5WjdVpVRmyPnJueIJvUCsmRZHq5/HJa9c5PXG+6LnZsnlprq6RDI/s7iY21esoOn009lxyimcU2J1VPp6XR1fq8vN969UrplIP/TDGGPOmsY4UiZRQ2f0mRaTxXu6HAwG8RXG+6LnXgv9Zw0NfGGftQj2eyor+c3q1TjsVvgct5P7l3n5yO4I93W5uG7fPhZ7HLy3en46Q1ZKjWHWtNDz7bnERuuDHjfSfC651kJ/pL2dj+20Lnu8pbSUO1etwiFCX98r7NhxOc88U8PLG1Zzedd5nMhGAK7Y8TJPvPp++vu3pDN0pdRRzKKEHp86t+yox8cTelckQsQdb6HXY0wshVHOnMZgkEu3bycGnFRYyP3HHks4sIPXXnsHGzaspbHxp4TDLQB4xMGXHN+jgF56KOLGjoW8+OLr2Lv3CzlzPpTKJZMuuWSLeELPMz0AuFxHH6B6+OCiWgCMCRMKNeH11qQoyplhjOGy7dtpC4cpdbm4b9VyWuuup67uJsAqTeXlraK6+qOUl5+P378UEQffqt/Ff+w5xKO8hXfyB6j7OgMD21i16rc4nb70vimlVELut9DtGnqeseYLG7449HDVHg/x/hxtZqjengs9XX7V1MT/dXYCcOviElq2vJG6uq8CUXy+pRx77P2cfPJm5s//FHl5yxPrrn503jLWFRQAcJ/3OgDa2h5k+/ZLtaWuVAbJ/YSeaKFbS8mNldDdDgdz3G4AmiIOXC6rRJPtk3S1hUJ8Zs8eAC4qjjF/35n0978COFmw4AuccsoWKisvOmLxbACHCF9ctAiAJ4Jz6a/9PgCtrfeyd+8XZuotKKXGMGsSuj9qLW4xVkKHYT1dfLkxje41e/fSHolQ5Ijwvu53E4l04HZXsnbtYyxZ8lUcDu9RH/+O8nLW5OcDcEfwbGprPwFAff036Ox8LOXxK6XGNnsSuhl/Qk/u6TK00EX2JvRnu7v5RVMTAJfFbqWMdgoLT2Lduo2UlJw5rudwiHDtAutcPNDWhm/e1yksXA/A9u0fIhweawp8pVSq5XRCj8Ri9MfnQre7LTqdY8/ae/ho0UUADA7uS02QKRY1hit37QJgOTu5gD9SWXkJa9c+ic83sX7lF1dWUul2EwV+3tTCqlV34nD4CQbr2b//+ukPXik1ITmd0HvsC6IwNNPiREsufv9SAAKBPSmIMPVua2jgpT7rvX+S71FT9T5WrfoNTqd/ws/ldTj4SHU1ALc3NuL2LWPRousBaGi4lYGBHdMWt1Jq4nI6oXcnzYU+NLBoAiWXpIQ+OLgPY6JHe1jGaQuF+MKerQC8jUc4Z+5prFz5y6NOTjaWy6urEaAxFOLh9nbmzfskPt9ijImwZ89npylypdRkzLKELrhcRWM+Lj6fS1MohMtrrVxkTJhg8GBK4kyVq7Y8QnfMTQG9XFfZxzHH3I6Ic0rPucjv562l1vTDv25uxuHwsmTJzQC0tz9ET8+LU45bKTU5uZ3Qk0ou+fTjdBaO2C1vuHjJJQZ0O2vB7pkeCOxNRZgp8cTBv/Kb7kIA/jPvOU5bdcu43vt4vL+qCoCH2tvpCoeprLyY/PzjAair+9q0vIZSauJyO6HH+6BLDCexcZVbYKjkAtAYBq/XWo4uW+rovb2vcPXubRgcLJYmrj/hGhwO97Q9/4UVFeQ5HISM4f62NkSEhQv/G7AGHPX1bZ6211JKjd+sSOjxudDHm9ALXS4KnVZp4vA6euYn9HC4g9tfvZ6NrAXgu8esw+ceu8w0EQUuFxdWWKNof91sLbxdWXkRfv8xANTXf2taX08pNT6zJKFb64qON6HD8MFF2dHTxZgor2x5P98PXwjA2YUuLqhalpLXeq9ddnmiq4vmUAgRJ/PnXw1AS8vdhEKtKXldpdTocjuhJ+ZCtxa3mEhCH6mnS6Yn9H37vshdXV7qWIgDw/eOWZuylYbeXFpKkdOJAf7Y1gZAVdV7cblKMCZIY+PPUvK6SqnR5XZCT6xWZC1uMdZMi8mSVy7y+5cAVkI3Gbp4cmvr/bxa9yPu4DIA/r26hjX2hFqp4HU4OL/cmlv+93ZCdzrzmTv3wwA0NPyIWCwy6uOVUtNvViT0ifRBjxup5BKNdhOJdExzlFPX1/ca27Z9kNv4GD0UU+x08pXFi1P+uhdVVgLwWGcnXWGrrFVb+3EAgsGDdHb+NeUxKKWGzIqEnmesxS0mVXIJhfD7h+rQAwO7pjHCqQuHO9i8+UJeji3lL5wHwNeXLKHK/gsjld5WVobP4SBsDA93WF90fv9SSkreBEBj4y9THoNSakhuJ3S7hp5PNzD5GrrLVYzHYw15HxjInCXYYrEIW7e+h97BOm7h0wCcXFjI5TUzsxBHvtPJ28qs6YV/3zp0EbS62iq7tLf/US+OKjWDcjuhJ1ro8cUtxl9Dj5dcArEYHZEI+fnHAmTMmprGGHbvvorOzke5j3exn4U4gJ+sWIEzRRdCR3KR3X3xkY4OBuwv0IqKi3A6izEmTHPz/85YLErNdrMioftj4586N26xb2hptb2BAPn5xwHQ3z+1QTORSC8HD/6QXbuuYufO/6Sp6deEw+0Tfp4DB26koeFWmqjiTvkIAFfW1nJiYeGU4puo88vLcYkQiMX4i112cTr9zJnzbwA0N/96RuNRajabJQnd6oUxkYRe6nZT5rImsdo9TQm9qelOnn9+Cbt3f4JDh75HQ8OP2L79Azz77Hx27/4UoVDzmM9hjGHfvuvZv/9LGODHnpsZNE5qPJ4ZuRA6XKnbzZtKrPOaXHapqvp/APT1bWRgYPeMx6XUbJTbCT1RQ594LxeAZX5ritndgQB5eVbJJRRqJByeeE+XQ4d+wvbtHyQcbsPh8FNW9i+Ulp6Lw5FHLBbg4MFbeO65pezbdz2RSO+IzxGJ9LJ9+4c4cOAGAF4q+BRPhqxFJ763bBlFrvSs+R3v7fJQezthe/754uIzEtcdWlt/l5a4lJptcjahR42hLzGwaPxzoSeLJ/Q9gQD5+asT2ydaR29tfYBdu64AoKTkbNav38vrXvcnjj/+L5x+ejNLltyM211BLNbPgQM38Pzzy9i//wb6+l4jGh1gcLCOQ4du5cUX19Dc/CsA/BWX8e3QxQD8S1kZF9tJNR0usPuj90SjPNPTA4CIk8rKdwPQ0nJP2mJTajbJ2YTeM8Jc6ONZrSjZ0qQWustVhNdrrS86kbJLKNTCzp2XA1ardc2ah/B65yb2u1wFLFjwWdav382CBV/A4fATDrewf//1bNjwOv75z3yee24hu3ZdSTB4ABEvy5bdwl2ez9MYCuN3OPjh8uUpGxE6HtVeLyfYg5j+3D50PWDOnPcA0N//Kv3929MSm1KzyZgJXUROEpEHROSAiAREpElE/iIip89EgJM18uIWE0voySUXYFI9XXbu/DjhcBsuVwmrV9+N05k/4nEuVzFLlnyV9et3MX/+5/B4ag/b73D4qa7+KCef/Cp7Cz7IjxoaALhh0SIW+Se++tB0e7vdffHPHUPlqKKiUxNrsra2aitdqVQbTwt9KeACfgpcCXwTmAM8KSJvSWFsUzJ8LnSHI3/CU8jGE3pzOExvJJK4MNrXt2lcj+/qepK2tvsJ46Kp+jaeGshjx8DAUacP8HprWbr0G5x2Wh2nnnqAE054ipNP3sYZZ3RxzDG3E3Qv5n3btmGAEwsKuGrevAm9p1R5u1122dzfT93gIAAiwpw5Q2WXTJ02QalcMeZVNGPMPcBhzSsR+TGwF/gk8LfUhDY1w1voLtfcoxw9sqVJLd+9g4PUFFmr3Pf2biAaHcTp9I32UIwx7Nnzee7jYu6WD9JeXwj1rwCwOi+Pz8yfz2Vz545aKhFx4PMtwOdbkNgWNYYPbt/OwWCQPIeD365ejduRGVWz9UVFlLpcdEYiPNLRwcfswU2Vle+hvv5bDAxso79/MwUFa9IcqVK5a1LZwBgzALQCE7vKOIMSXRYlhovohMstAHPcbgrsedF3BwIUF78BAGOC9Pa+cNTHtrf/iVt7l3IrV9JuChHi6x7B1oEBPrxjBx/avj0xGGc8Pr93Lw/ZNepbly9nRV7ehN9TqjhFODdedkmqoxcWrsPnsyY304ujSqXWuBO6iBSKSIWIHCMiXwOOAx4b5diuo92AiWfXCZrKXOhxInJYHd3jqSQvbxVglVNGY0yU7+z8E7/EGgJ/UUUFB049lciZZ/LiiSfyr/boyl81N3PChg0829191DgGo1E+tmMH36yvB+CqefO4rLp6wu8n1eJ19Ec7Owna3Retsot1cbS19d60xabUbDCRFvovsVrl24GrgZ8AGbuAZGIudAkBk0voMFRH3zUwAEBx8Rut5+8ePaE/VXcv3w5dBMDbil3cvXo1830+HCKcVFTE/ccey7eXLsUtws5AgNNffpkPbtvGAbv2nGxTby/rX3qJ2xsbAbi4ooJvLV06qfeSaueWlSHAQCzGk11die0VFda5CAR2am8XpVJoIiNRbgBuA+YBHwC8gBvs1SOSGGOOmj1nopU+NBe6lSQnm9CPzcvjPmBTn9WXvaTkTBobb6O7+xlisfARF1pD0SCX7x8khJdqRy93rznviDq3iPDp+fM5p7SUy7ZvZ1NfH3c2N/OblhbOKytjgdeLQ4TX+vt5squLGNY371cWL+baBQtwpLGL4tHM8Xg4ubCQF3p7+XNHB2+xW+yFhevweGoJhQ7R3v4H8vNXpjlSpXLTuFvoxpjXjDF/M8b8EjgXWAfckarApiqR0IkvbjG5hL7Onhvltf5+QrFYoo4ei/XT27vhiON/svNetptFCDHuWLGY4qOM3jy+oIAN69Zx+4oV1Ho8RIzhofZ2bm1o4AeHDvEPO5mv8Pt55sQT+cLChRmbzOPivV2S6+giQkXFOwFrEWmlVGpM9qJoGPgDcJGIpL8T9AiGFreIjxKd3B8E8YQeMobN/f34fPMSdfSmpl8ddmwo0sd3mq2Ee66vnrfOPW7M53eK8NGaGvaceio/P+YY/qOmhgsrKvjXigour67moeOO47WTT2Z90fQu9Jwq8Tr6zkCA3XaZCkgk9J6e5wkGm9ISm1K5bip93vxYHTdmdnq/cYrX0POMNRR9si30Gq+XufZiERt7rTlWamr+A4Dm5rsIh4dqxbfvuJMDWAOC/mf5GRN6Ha/DwYerq/nxihU8cNxx/P6447jtmGM4v6ICT4Z0TRyPdYWFVLqtMtQjSYOMSkrOwuksAgzt7Q+lKTqlctt4RooeMUmIiBQBlwD1xpiWVAQ2VUNzoU98cYvh1tnD2uMJfe7cD+Jw5BOLDdDUdAcAnV3P8p1W64+Vt/qbOLk8My9cpppDhPNGGDXqcHgoL387oGUXpVJlPE2/e0TkzyJynYh8RERuADYDC4DPpDa8yUuUXEwnMMWEbpddNvYNlW/mzv0AAAcO3EB9/Xf4yeZvsg9r+tobV54z6dfKBfE6+uOdnYf1sy8vt8ounZ2PEYn0pSU2pXLZeBL6r4E84L+AHwMfB14BzjbGZOy8qEcubjH5TjXxhP5qXx8hu3/1/Pmfw+2eQyTSxe49V/OLiLWe53nFHk4uTt/Mh5ngraWlOICgMTye1H2xvPw8RNwYE9QFpJVKgTETujHmF8aYs4wxc4wxbmNMpTHmHcaYJ2YiwMlK1NCZ+ALRwyVfGH3JLrv4/YtZt24jhYXreZbT2M1yAK5fOvaF0FxX6nZzerH1BfpIUm8Xl6uYkpKzAC27KJUK2XO1bYKGui1Obi70ZLVeL8faw+zvTVqVx+ebxwknPM0DBd8F4NzSUk7Jkt4oqRavoz/c0XHYpFwVFRcC0N7+MLFYZMTHKqUmJycTetQYeqe4WtFw/zZnDgD3tLQQS0pQ/9fZxYY+q6/7lxYtmtJr5JJ498X9g4PsSOq+WF7+DgAikU56ep5OS2xK5aqcTOi907C4xXDvsRP6oVCIp+y5V4KxGNfu3QvAm0pKEmUGZQ2aqra7eyb3dvH55lNQcAIAbW3afVGp6ZSTCX34XOgi3qNOdTsey/PyEt0Xf2rPq3Lt3r280t+PA7gxDQs0ZzIRGVr0IqmODkOt9Pb2P+gc6UpNo9xM6Ekt9AL6plxuifvgXGtO9V83N3PWyy9zy8GDAHxx4UJO09b5EeLdF5/s7j7sr6aKigsACAR2MzCwIy2xKZWLcj6hW4tbTE9Cv6KmhnfZizE/YZdd3lxaynULF07L8+eaN5eW4hIhbAyPdXYmthcUnIjHYy2AoaNGlZo+OZ3QfRK1F7eYnoTucjj431WreO+cOSz1+fjx8uU8smYNriwamj+Tilwu3mD/5ZJcRxeRpLLLH9MSm1K5aCLT52aNeA29UMJgpjaoaDiPw8H/rl49bc+X695eVsbjXV38ub0dY0xiyb2KigsS0xCHQm14PBVpjlSp7JeTTcuhudCntriFmrp4Hf1QKMRr/f2J7SUlb8LhyANidHT8OU3RKZVbcjyhT20udDV1q/LyWOj1Aof3dnE6fZSVvRWAtjYtuyg1HXI6oedjDWjRhJ4+IjK06EVSHR2gvNzq7dLZ+VdisSMWvlJKTVBuJvTEKNGpz+Oipi7eH/2Z7m46w+HE9vLyfwGEaLSPrq5/pCc4pXJIbib0xFzo8cUttI94Op1dWopHhCjwt6Tuix7PHIqKTgV01KhS0yHHE7o1dau20NMr3+nkrBLr/+DIUaNW2aW9/Y86alSpKcrxhD71xS3U9IjX0R/p6DhscrOKCqs/ejBYT1/fK2mJTalckZsJPbGe6NSXn1PTI15HbwmHeblvaLWivLzV+HxLAB01qtRU5WZCP2IudK2hp9vyvDyW+a01V5PLLjpqVKnpk9MJfbrmQlfT4+0jLB4NQ5N19fZuIBhsmPG4lMoVOZfQYylY3EJNj3gd/fmeHtpCocT24uI3JOarb2//U1piUyoX5FxC741GiV9ys+ZCd9lDzFW6nVlcjN/hwAB/Teq+6HC4KS+3FtnWOrpSk5dzCX2kudDjE0Kp9PI5nZxTWgqM3n2xs/NRotH+Ix6rlBpbTif0fPqnvPScml7xOvojHR2EY7HE9rKytwFOYrFBOjsfTVN0SmW3nE/oWj/PLO+w6+idkQiPd3UltrvdpZSUvBHQUaNKTVbuJXT7gqhXoriJaELPMPN8Pk4vKgLgvtbWw/YNjRr9E8bEjnisUuroci+h2y30QrEmgdKEnnniy/g90NZGJKnsEh81Gg4309v7YlpiUyqb5WxCL5BBQAcVZXPTTPUAAB5oSURBVKKL7YTeFg4n1mYF8PuXkpdnrQalc6QrNXFjJnQROVlEbhWRrSLSLyJ1InK3iCybiQAnamhQkS5ukakW+HysLywE4N6WlsP2xQcZtbX9fsbjUirbjaeFfg1wEfAo8EngduAs4GURWZW60CYnXkMv0EFFGS1edvl9WxvRwybruhiAgYHt9PdvSUtsSmWr8ST07wALjTH/ZYz5mTHmRuANgBsr2WeUxEyLurhFRosn9NZwmH8m9XYpLFyHz7cIgJaWe9MRmlJZa8yEbox5xhgTGrZtF7AFyLwWerzkojMtZrRFfj8nxcsuSb1dRITKyncB0Np6X1piUypbTeqiqFhDL6uAtlH2dx3tBqTsSmU8oftj8cUt9KJoprpklLJLPKEPDGyhv39bWmJTKhtNtpfL+4Ba4HfTGMu0GFpPVBe3yHTxsktTKMTTSb1dCgtPwetdAEBrq5ZdlBqvCSd0EVkJ3Ao8Bdw10jHGmJKj3YDukR43HXTq3OyxxO/nxIIC4PBBRlp2UWpyJpTQRWQu8DDQCVxiMnA4nyb07BJvpd/f2nrY0nTxhN7f/xoDAzvSEptS2WbcCV1EioFHsOrf5xpjmlIW1RR06WpFWSWe0BtCIZ7t6UlsLypaj9c7D9DeLkqN17gSuoj4gIeAFcD5xpiMbDJFYjF67Bp6Ib2A4HQWpjcodVTL8/I4Pj8fgHuSBhmJOBJ90ltbM+5SjVIZaTwjRZ3APcBpWGWW51Ie1SR1DZsL3eksQiTnZjfIOf+vqgqA3zQ3E0qa22XOnH8DrLJLX9+raYlNqWwynmz3beACrHJLmYi8P+l2YWrDm5jOpIReRA9ud1kao1Hj9YGqKhxAeyTCQ0kLXxQVrcfvXw5Ac/OI19+VUknGk9DX2vfvwOrVkny7JUVxTUpyQi+kF5erNI3RqPGq8Xp5m73wxS8bGxPbRYSqqvcD0Nz8G4yJpiU+pbLFeEaKnmWMkVFui2YgxnGLJ3QXMXwMags9i3xo7lzAWsmoIRhMbI8n9FCogc7Ov6clNqWyRU4VmDvD1hzoRY4QArhcmtCzxTsqKih3uYgBdzU3J7b7/UsoLj4D0LKLUmPJqYTeYbfQiyQ+da6WXLKF1+HgffbF0V82NmKS+qRXVX0AgNbW3+sC0kodRU4l9HjJpdAeVKQll+wSL7vsCAQO65NeWXkJIh5isX5aWx9IV3hKZbzcSuh2yaUgMXWuttCzydrCQk6wpwL4edLFUbe7lPJya3k6LbsoNbrcSujxUaL21LnaQs8+/15dDcBvW1pot7+gAebOvRSAzs6/MTh4IC2xKZXpcjShW32Z9aJo9rm0qooip5NALMbPklrpZWXn4fHUAIaGhtvTF6BSGSzHE7qWXLJNocvFh+1W+q2HDhGxR446HG6qqz8CQGPjz4nFQqM+h1KzVW4l9EQN3ZqYS0su2enK2loEqA8GebBtaA2V6uqPAg7C4Wba2h5MW3xKZaqcSugdiV4uelE0my31+zm/vByA7x86lNju881LXBxtaPhJWmJTKpPlVELvHJbQtYWevf6rthaAf3Z383Jvb2J7be0VAHR1PU5///a0xKZUpsqZhB6Oxeizp84toA8RDw5HXpqjUpN1Tmkpq/Os/7/vHTyY2F5a+hZ8viUANDbelpbYlMpUOZPQu0aYmMtay1plIxHhk/OsBS7+t6WF/YGAvd1BTc3HAGhqukNHjiqVJGcSuk6dm3suraqixuMhYgw31dUlts+d+yEcDh+RSBeNjT9PY4RKZZacTOgF9Gkf9Bzgczq5dsECAH7R1ETd4CAAHk8lc+deBkB9/beJxcKjPYVSs0rOJPQOu8uiOzF1rvZwyQUfra6m2uMhPKyVPm/e1YCDYLBOl6hTypYzCT3eQi9yBHXq3Bziczq5xm6l/7yxkYN2Kz0vbxmVldaao3V1Nx82O6NSs1XOJfRCdOrcXHN5dTVVbjehYa30+fM/B0B//6t0dPw1XeEplTFyJqHHSy5FoqNEc43f6eRzdiv9p42N7LF7vBQVnURJyZsAqK+/OW3xKZUpciaht9gJvYROQEsuueaKmhrme72EjOEze/Ykti9YcA1gDTTq7n4mXeEplRFyJ6GHrMmaSow194deFM0tfqeTm5dYA4oebGvjsU7ri7u09C0UFp4EwN69n9dauprVciaht8ZLLjFrPUq3uzKd4agUeM+cOby+qAiAj+/cyWA0ioiwePHXAejuflJr6WpWy5mEHm+hF9MBgNs9J53hqBQQEX60YgUuEXYGAnzNvkBaVvbmRC19374vYEwsnWEqlTa5k9DtFnqpXUP3eLSFnoteV1DAZ+bPB+Cmujq29ltD/5cssVrpfX0v09p6b9riUyqdciKhR41JLFdWQhegJZdc9qWFC1nq8xE2hst37CBmDEVFp1BR8a8A7Nv3RR09qmalnEjoHeEw8T+yS+jC6SzG4fCkNSaVOn6nk5+sWAHA0z093NbQAMDixTcCDgKBXTQ0/DiNESqVHjmR0FuSFhMupROPR+vnue7NZWVcWlUFwNV79rClv5/8/NWJZer27fsiwWBTOkNUasaNK6GLSLWI3CQij4tIr4gYETkrxbGNW2toaH3JYrr1gugs8b1ly1jk8xGIxXj3li0MRKMsWfI1XK4yotEe9u79XLpDVGpGjbeFfgxwDTAPeDV14UxOYlCRBHES0wuis0SJ283dq1fjEmHrwABX7d6N213OkiU3AdDcfBddXU+mOUqlZs54E/pGoMIYsxz4ZgrjmZR4l8VSxwCgXRZnk/VFRXxt8WLAmhbgt83NVFf/O4WF6wHYtes/9QKpmjXGldCNMb3GmPZUBzNZrYkuiz2A9nCZba6eP59zS62Rwf++YwcbevtYseJHgNDfv5n6+m+lN0ClZkhKLoqKSNfRbkDxdL5evOQSH1SkF0VnF4cIv161iiV2Pf2CzZvpdK+mtvZKAPbv/zJ9fRlXKVRq2uVGL5f4KFHTAmgLfTaq8Hh4eM0aSlwumkIhzn/tNSoW3IjfvwxjwmzbdimxWGjsJ1Iqi6UkoRtjSo52A7qn8/XiJZfimNVNTVvos9PK/HzuP/ZYXCK81t/P/9u+j8UrfgU46O9/hf37b0h3iEqlVG610BOjRDWhz1ZvKi3lNnvQ0V87O7niUBE1868FoK7uJrq7n01neEqlVG4k9GHzuGjJZXb7cHU1X7d7vvy+rY1P978bR95JQIytW/+NcDhjr+8rNSVZn9DDsVhi+bmheVwq0hmSygDXLlzIFxcuBODhjk4+xXfpklqCwTq2bXu/zsioclLWJ/SDwWDi50pacbnKcDhcaYxIZYr/WbyYW5cvxwG8NBDhk647OMACOjr+woEDN6Y7PKWm3bgTuohcJyLXAZfYmz5gb7syNaGNzwF7FXiwErpeEFXJPl5bywPHHYff4aAu7OI/5Wc8xevZv/96XQxD5ZyJtNC/Yt/ea//7w/a/PzPdQU3EAbuFPscZxEMYt7sqneGoDHRBRQVPrF1LjcdDv3HzRW7k+1zJS5s/QH//1nSHp9S0GXdCN8bIKLdFKYxvTPEWerXD6gnp881PZzgqQ51cVMTGdet4Y7E1pu0BLuIjsZu575UrCIVa0hydUtMj62vo8YRehfVL6fVqQlcjm+v18ve1a/nq4sW4gAMs4t9DX+SyDT+ldbA33eEpNWU5k9DnxOoBTejq6JwifGHhQp498USWeqJEcfHb0OtZ/vzTfL++jlBMe7+o7JX9Cd2uoVdE9wDg8y1IZzgqS5xUVMSWU8/m82WHyKOfbuPjk3v2svKFF7i9oYH+aDTdISo1YVmd0GPGUBdvoXMI0Ba6Gj+vw8FX17yXv9c8x/k8hIMo+wYH+djOndQ88wyf2LWLTb29GGPSHapS45LVCb05FCJk/7LNxZrHRRO6mggR4ZTl/8O3awa4g8t4K3/FQ5SeaJQfHjrECRs3svz55/ncnj08191NTJO7ymBZPQInuQ96Fc04HPm4XCVpjEhlIxFh+fJbiUR6mN9yEx/nRzxf8k3uCx7PrkCAPYODfLO+nm/W11Pj8fCvFRW8s6KC1xcXk+d0pjt8pRKyO6Hb9fNiR4S8WACfbyUikuaoVDYScbBy5R0YE4LW+3hr18f48LzP0lf1RR5oa+P+1la2DAzQEApxa0MDtzY04BJhXUEBZxQX84aSEs4oLqbc7U73W1GzWFYn9Hj9vMbZDzEtt6ipcTjcrFr1W0Q8tLT8hoMHv0mtCXL9su9yw+LF7BgY4IHWVu5va2NDby8RY3i+t5fne3v59sGDAKzOy+MNxcWJJL/A69VGhpoxWZ3Q98YHFYm1UpEmdDVVDoeLVavuRMRNc/OvOHTo+4TDbaxc+UuOycvj2oULuXbhQjrCYZ7p7uaf9m1Dby9hY9g6MMDWgQFua2wEoNDpZIXfz4q8PE4uLOT04mJW+P2UaktepUBWJ/SXe63BIIs5AGhCV9NDxMnKlb/A6cynoeFHtLT8hnC4hWOP/T0uVyEAZW4351dUcH6FNbNnIBrlhd5e/tnVxVPd3TzT00NvNEpvNMrGvj429vXx25ahEallLhdL/X6W+f2J+2Pz8lidn691eTVpWZvQI7EYr/T3A7AsZq0XqcP+1XQRcbB8+Q/xemvYt+86OjsfZdOms1iz5o94vbVHHO93OjmzpIQzS6yL8pFYjG0DA2wfGGBnIMDW/n6e6u6mzr7u0xGJ0NHby4u9h49QFWCJz8fyvDyW+Hws9vlYkZfHqrw8lvj9OLV8o44iaxP61oEBBu1RfUsizwHaQlfTS0RYuPC/8Xiq2bHjcvr6XmLjxlM47rg/UFR00lEf63I4WFNQwJqCgsO2d4TD7AkE2B0IDN0PDrJjYIDWcBgD7BkcZE9SD644jwjL7fJN/D5ezpnjdmutXmVvQt9ot2wqXA4qI1a90u9fls6QVI6qrv4wXm8tW7a8m1CogU2b3sjKlXcwZ867J/xcZW43ZW43JxcVHbGvJRRic38/W/r72Ts4yD472e8aGCBoDCFj2DIwwJaBgSMem1yrX+H3szx+7/dTovX6WSPrE/oa3yDSBw5HHj7fovQGpXJWWdm5nHjic7z22vkMDu5l69b30Nu7kcWLvzptC6rM8Xh4k8fDm0pLD9seNYZ9gQDb7PLNTvt+18AAh+z1dJNr9cNVuN1Wrd7no8brZb7Xm0j8C3w+LePkkOxN6PYHd5XTGiGan78akawe+KoyXH7+Ktate4EtW95FV9c/qK+/mZ6e51m9+m683rkpe12nCMvy8liWl3fEvr5IhN2BgJXgk5L9zoEBOuylGdvCYdrCYZ7r6Tni8V4RliWVcZYl3Wq9Xhya7LNKVib0SCzGK3ZCX262AZCXd2w6Q1KzhNtdzute9zf27ftv6utvprv7CTZuPIGVK++irOzNMx5PgcvF2sJC1hYWHrGvPRxm18AAewYH2R0IsDcQoDEU4sDgIHsHB4kYQ/AoZRyvCAt9PuZ7vczzepnv8zEv/rN9X+pyae0+g2RlQn+yu5uAfUF0cegpAPLzNaGrmeFwuFi69BsUFZ3G9u2XEQo18eqrb2HevE+zZMnXcDi86Q4RgHK3m/LiYk61F/VIFo7F2Dc4eFiLfqd9obY+GMQAQWOsfYHAqK+R53AckeSrvV7mejxUezzMtW/aFXNmZGVC/8Eha2bF04sKKel9DgPk5x+X3qDUrFNZeSEFBS+xbdv76Ol5joMHv0Nn56OsWnUnBQXHpzu8o3I7HFYdfYQyzmDUmnVydyDAgcFBDgaDHAwGqbfvDwaDiUnxBmKxMZM+QL7DQZXHQ5XHwxy3+4if53g8VNk/l2irf9IkHVODikhXcXFxcVdX14Qfuy8QYOnzz2OAXy0tZMEeq/vYqace0LnQVVrEYhEOHLiRAwe+AsQQcTF//jUsXHgdTqcv3eFNu5gxtIXDRyT5ejv5N4VCNIVCdE9yTnm3yBGJfo79BZD8c4XbTYnLRZ7DMWu+AEpKSuju7u42xow4C2HWJfSrd+/mOwcPMs/r5bkl9eza9h6czkLOOKN71vynqszU3f0s27d/iEBgBwB5eStZvvzHlJaeld7A0iQQjdIUCtEcCtEcDtNylJ/jF3Anwy1CictFqctFictFhdvNQp+PCrebIqeTQpeLIqeToqT7UvtW4HRmVd4YK6FnVcklZgx/am8H4OM1NYQG/gJY9fNs+k9Ruam4+DROOmkTBw58hfr6mxkY2M4rr5xNZeW7WLr0W/h8C9Md4ozyO50s9vtZ7PePeWwoFqM1KdG3hMNWwh/h55ZQiOS2f9gYWsNhWsPhCcfoBOvLwG7tl7pcFLtcFDud1v2wfxeNsM/ryJzedVnXQh+MRvldaytvLyujfstb6e5+gpqaK1ix4kcpiFSpyenre4WdO6+gp+dZABwOH/PmfYr58z+L2106xqPV0cSMoSMcpiMSoTMSoSsSoTMctu4jEZpDIQ4Eg3SGw/TY8+n0RCL0RqOJzhTTySsyarIvdjrxO50YYyh1u5mbdK1gbUHBhBuiOVdyiYtEenj66XKMiXDccX+gouKCaYxQqakzxtDS8hv27PkcoVADAE5nMfPnX828eVclJvpSMycci9ETjdIZDg99GQz7YuiORumORIZu9r97IhF6olGmI2MWOp30vOENE35cTpVcknV2PoYxEUTclJScne5wlDqCiFBV9T7Ky9/JwYPfpr7+20Sj3ezf/yUOHvwONTX/QW3tJ/B6a9Id6qzhdjgodzgmvRBJzBj6khO+3fpPTvzD9w3GYhiscQHx6wZVHs/0vjFb1rbQd+z4GI2Nt1NScjZr1/59GqNTKjXC4Xbq6m7m0KEfEItZ3fxE3MyZ8x6qqz9KcfEb9FrQLBD/UihyTbw9PVYLPXOq+RNgjKGj4xEAysrOS3M0So2P213O0qXf4NRTD7Bo0Q243ZUYE6a5+dds2nQmL7ywkgMHbiIQ2JvuUFUKOUQmlczH9dzjOUhEvCLyDRFpEJGAiDwnIuekJKJx6Op6gmCwHoCysrelKwylJsXjqWTRoi9x6qkHWLHipxQWngJAILCTffs+z/PPL2XDhnUcOPA1entfwpjpv5CnctO4Si4i8lvgYuAWYDdwGXAScKYx5tkJv+gUSi7GGF566VR6e1+gqOh0TjjhKf0zVWW9vr5XaWz8Oa2tvyMUajpsn9tdQUnJORQXn0Zh4XoKCtbm5IAlNbYp93IRkVOA54FPGWNusbf5gM1AgzHmjRMNaioJvaXlPrZuvQSAtWv/SUnJGRN+DqUylTFRurufobX1XtrbH2Zw8Mjyi4ibgoLjKSg4gby8Y/D7V5CXdww+32IcDp37PJdNR0K/GbgKKDPG9CVt/zzwVaDWGNM4kaAmm9BjsQgvvriaQGAX5eUXsGbNHyb0eKWyTSCwh46Ov9LV9QS9vS8wOLj/KEc78Xpr8Hhq8Hpr7Z+rcblKRrgV43D4cTh8OBxeRHTyrGwwHQn9b0CVMeZ1w7afAzwKvN0Y88iwfWNl6uLi4mImmtCNMXR2/o29e7/AqlW/0hkW1awTCrXQ0/MCvb3P09+/lYGBHQQCuzAmNKXnFXHbyT2e4K0kb60xYN1bSf/w++T9o12SG70kOtL28R87Pc+bHi5XCatX/++EHzcd/dCrgUMjbI+3ymesE62IUFb2VkpL36J1czUreTxzqKg4n4qK8xPbjIkyOFhHILCTYPAQwWADoVCDfd9EJNKVuBkTHPF5jQkTjYaJRntH3K+ml9tdlZLnHU9C9wMjfQoGk/YfZrRvjzi7BX/kJM3jpMlcqSEiTvz+xfj9i8c8NhodJBrtJhLpIhYbHHYLHvazMVEgdti91eNmpG1RRv5rf7QKwPiPTdXzppPTmZpRwuNJ6AFgpBn7fUn7lVJZwOn04XT68HhS00JU6TWefuiNWGWX4eLbGqYvHKWUUpM1noS+CVgpIgXDtq+371+Z3pCUUkpNxngS+n2AG/hIfIOIeIEPAU8bY7SFrpRSGWDMGrox5nkRuRe4WUSqgT3AB4GFWCNGlVJKZYDxzhBzKfAV+74UeBWr//nTqQpMKaXUxIwroRtjBoHP2jellFIZKF3zoccAKS6edFd0pZSadbq7uwGMMWbE65/pSugRrAuyPaMcEs/03TMTUc7Q8zY5et4mR8/b5EzlvBUBMWPMiNWVtCT0scTnghlrxKk6nJ63ydHzNjl63iYnlectK1csUkopdSRN6EoplSM0oSulVI7QhK6UUjlCE7pSSuUITehKKZUjNKErpVSOyMh+6EoppSZOW+hKKZUjNKErpVSO0ISulFI5QhO6UkrliIxJ6CLiFZFviEiDiARE5DkROSfdcWUKETlLRMwot5XDjj1dRJ4SkQERaRKR74lIXrpin0kiUi0iN4nI4yLSa5+fs0Y59gIReUlEBkWkTkS+LCJHzGInIiUicruItIpIv4j8XUTWpvzNzKDxnjcR2T/KZ/CmEY7N6fMmIieLyK0istV+f3UicreILBvh2HH9Tk41D453xaKZcAdwMXALsBtrebtHRORMY8yzaYwr09wCbBy2LbGuq/0L8xiwBfg0MA/4DLAEeMcMxZhOxwDXYH2GXgVOH+kgETkPeBD4O/AJYA3wJaDC/nf8OAfwsL3/W0A78HHgHyKyzhizJ2XvZGaN67zZNmJ9DpNtTv7HLDlv1wCvB+7FOmdzgSuBl0XkFGPMNpjw7+QdTCUPGmPSfgNOAQxwVdI2n/2Gnkx3fJlwA86yz9GFYxz3Z+AgUJC07SP2Y9+U7vcxA+epECi3f77Qft9njXDcFqzE5EzadiMQBZYnbXv38PMOVAKdwJ3pfr9pOG/7gQfH8Xw5f96wvvQ8w7YtBwaBO5K2jet3cjryYKaUXN4FhIGfxTcYa9m7nwNn2ItTK5uIFI5SGigC3oL1C9OXtOtOoA/rlyynGWN6jTHtRztGRFYDq4HbjDHRpF0/wipDXpy07V1YfwH9Iek1WoHfAReKiHu6Yk+n8Zy3ZHZp4GhlvJw/b8aYZ4wxoWHbdmE1FlbBhH8np5wHMyWhnwBsH/aGAV4ABMiZuts0uAtrpaeAiPyfiKxJ2rcGq4y2IfkB9oduE9Z5VkPnYfh5asBqSZ0w7NiNxm4uJXkBq1V7RL10Fngr0A/0i8geEbl8hGNm5XkTEQGqgDZ700R+J6ecBzMloVcDjSNsj2+rmcFYMlUIuA/4JPBO4AasP9GeEpEV9jHxb/DRzqWeR8tEzpN+Ng/3KvBlrL9iPoqVuG4TkWuHHTdbz9v7gFqsv0Rghj9rmXJR1A8ER9g+mLR/VjPGPAM8k7TpjyLyENY3/5exPkjx8zTauZz159E21nnKG3asfjZtxpgLkv8tIr8EngK+KCI/NsbE18mcdefN7m12K9b5uMvePJHfySmfs0xpoQcA7wjbfUn71TDGmFeAR4F4t6b4eRrtXOp5tEzkPOln8yjsaxC3YH0Jnpa0a1adNxGZi9WrpxO4xBgTs3fN6GctUxJ6I0N/miSLb2sYYZ+y1ANl9s/xP81GO5d6Hi0TOU/62RxbvX1flrRt1pw3ESkGHgGKgXONMU1Ju2f0s5YpCX0TsFJECoZtX2/fvzLD8WSTJUCr/fNmIAKclHyAiHiwLqhsmtnQMlb8PAw/TzVYfYQ3DTt2nX2xK9l6rF4Ku1MVZBZZYt+3Jm2bFedNRHzAQ8AK4HxjzI5hh0zkd3LKeTBTEvp9gBurbyZgdYsCPgQ8bfc+mNVEpHKEbWcAZwN/BbDrl48CHxj2ofgAUIA1AGLWM8ZsAbYDl4uIM2nXFUAMuD9p231YF6PeGd8gIhXAJcAfjDHh1EecGUSkzB4wlLzNB3wW6AWSB77k/HmzPzv3YJWaLjHGPDf8mAn+Tk45D2bMfOgi8jusAQ3fBfYAHwROBs42xjydztgygYj8HRjAujDaBhwHXA50AycbY+rs4060j9mM1Z91HnA18Lgx5u1pCH3Gich19o+rgPcCvwD2AV3GmB/ax5wP/BFrpOg9WOfzSqy+6R9Pei4n1kWuY7FGPLZhjXicD6wzxuRESxPGPm8ichnw31iJZz9QjvV7ugK4whjzk6TnyvnzJiK3YPU6e4ihXi1xfcaYB+3jxv07OeU8mO7RVsNGRH0Tq440iNX38s3pjitTbsB/Ac9jDaEOA4ewfuEWjHDsGcDTWBdRmoHvA/npfg8zeK7MKLf9w467EHjZ/rzVY3UFdY3wfKX2L2IbVv/rx4ET0/0+Z/q8Aevs5HUQqzdGD/APrFLDSM+X0+fNfu/j/ayN63dyqnkwY1roSimlpiZTauhKKaWmSBO6UkrlCE3oSimVIzShK6VUjtCErpRSOUITulJK5QhN6EoplSM0oSulVI7QhK6UUjlCE7pSSuWI/w/LWBrHE2GT1AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "G = np.zeros(16)\n",
    "G[0:6] = [3.7e8, 1.5e8, 1.1e8, 1.1e8, 2.5e7, 1e7]\n",
    "lambda_ = 17\n",
    "rho = 1.6\n",
    "stacked = generative_model(G, lambda_, rho)\n",
    "plt.plot(coverage_range, stacked, \"y-\")\n",
    "plt.plot(coverage_range, P, \"c-\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6431330311733933"
      ]
     },
     "execution_count": 219,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func(17, 1.5, G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimizer\n",
    "\n",
    "Optimize one variable at a time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of optimizing all $G_i$, we bin the larger values together, in geometrically sized bins. Within the bin, the value is presumed to be the same. A total of 22 bins.\n",
    "\n",
    "- $G_1$ to $G_8$ (single copy to 8-copies) as usual\n",
    "- $G_{9}$ to $G_{11}$, size 3\n",
    "- $G_{12}$ to $G_{16}$, size 5\n",
    "- $G_{17}$ to $G_{23}$, size 7\n",
    "- $G_{24}$ to $G_{32}$, size 9\n",
    "- $G_{33}$ to $G_{45}$, size 13\n",
    "- $G_{46}$ to $G_{64}$, size 19\n",
    "- $G_{65}$ to $G_{91}$, size 27\n",
    "- $G_{92}$ to $G_{128}$, size 37\n",
    "- $G_{129}$ to $G_{181}$, size 53\n",
    "- $G_{182}$ to $G_{256}$, size 75\n",
    "- $G_{257}$ to $G_{362}$, size 106\n",
    "- $G_{363}$ to $G_{512}$, size 150\n",
    "- $G_{513}$ to $G_{724}$, size 212\n",
    "- $G_{725}$ to $G_{1024}$, size 300"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 385,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 11), (12, 16), (17, 23), (24, 32), (33, 45), (46, 64), (65, 91), (92, 128), (129, 181), (182, 256), (257, 362), (363, 512), (513, 724), (725, 1024)]\n"
     ]
    }
   ],
   "source": [
    "bins = []\n",
    "for i in range(1, 9):\n",
    "    bins.append((i, i))\n",
    "for i in (8, 16, 32, 64, 128, 256, 512):\n",
    "    a, b = i + 1, int(round(i * 2 ** 0.5))\n",
    "    bins.append((a, b))\n",
    "    a, b = b + 1, i * 2\n",
    "    bins.append((a, b))\n",
    "\n",
    "print(bins)\n",
    "# for a, b in bins:\n",
    "#     print(\"- $G_{{{}}}$ to $G_{{{}}}$, size {}\".format(a, b, b - a + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 470,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average based on last 5% of histogram [last_N=500]: 36\n",
      "[         0 1281576854   89292133 ...         36         21     322132] 10002\n",
      "[         0 1281576854   89292133 ...         36         36         36] 16704\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "# NOT USED\n",
    "\n",
    "# hist often end with the last bin having big numbers,\n",
    "# we fix this by taking all the kmers from last bin and\n",
    "# distribute it in all the smaller bins\n",
    "def process_hist_redistribute(hist, inflate=math.e):\n",
    "    N = len(hist)\n",
    "    last_kmers = hist[-1] * (N - 1)\n",
    "    newhist = [0] * (N - 1)\n",
    "    bonus = last_kmers / (N - 2) * inflate  # except 0 and N - 1\n",
    "    for i, h in enumerate(hist[:-1]):\n",
    "        if i == 0:\n",
    "            continue\n",
    "        newhist[i] = hist[i] + int(round(bonus / i))\n",
    "    return np.array(newhist)\n",
    "\n",
    "# hist often end with the last bin having big numbers,\n",
    "# we fix this by taking all the kmers from last bin and\n",
    "# extend the range\n",
    "def process_hist(hist, last_percentile=5):\n",
    "    N = len(hist)\n",
    "    last_kmers = hist[-1] * (N - 1)\n",
    "    last_N = N * last_percentile // 100\n",
    "    average_last_N = int(math.ceil(np.sum(hist[-last_N-1:-1]) / last_N))\n",
    "    print(\"Average based on last 5% of histogram [last_N={}]: {}\".format(last_N, average_last_N))\n",
    "    new_hist = list(hist[:-1])\n",
    "    i = N - 1\n",
    "    while last_kmers > 0:\n",
    "        propose = min(average_last_N * i, last_kmers)\n",
    "        new_hist.append(average_last_N)\n",
    "        last_kmers -= propose\n",
    "        i += 1\n",
    "    return np.array(new_hist)\n",
    "        \n",
    "new_hist = process_hist(hist)\n",
    "print(hist, len(hist))\n",
    "print(new_hist, len(new_hist))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 475,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17 1.6 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "Iteration 0\n",
      "16.535391641021622 1.824782682224368 [4.06580616e+08 1.90629217e+08 1.61240653e+08 9.98586715e+07\n",
      " 2.60892013e+07 1.39478571e+07 1.02173244e+07 7.22086328e+06\n",
      " 3.07962545e+06 1.25434492e+06 5.04331206e+05 2.34760837e+05\n",
      " 1.35340907e+05 6.00976883e+04 2.65464502e+04 1.27983140e+04\n",
      " 6.66397886e+03 3.39601283e+03 2.21022146e+03 1.18324862e+03\n",
      " 0.00000000e+00 9.99999970e+09] 4.763132413865311e+17\n",
      "Iteration 1\n",
      "Iteration 2\n",
      "Iteration 3\n",
      "Iteration 4\n",
      "Iteration 5\n",
      "Iteration 6\n",
      "Iteration 7\n",
      "Iteration 8\n",
      "Iteration 9\n",
      "Iteration 10\n",
      "16.469555084106965 1.4938343161324967 [3.63614539e+08 1.41661365e+08 1.03021640e+08 1.13733508e+08\n",
      " 2.12133057e+07 1.33481465e+07 5.99421404e+06 7.17084036e+06\n",
      " 3.04166688e+06 1.21734997e+06 4.90230895e+05 2.22096140e+05\n",
      " 1.30195533e+05 6.08372685e+04 2.63639725e+04 1.27764209e+04\n",
      " 6.66882778e+03 3.31330088e+03 2.03694922e+03 1.29709294e+03\n",
      " 0.00000000e+00 9.99999970e+09] 3372693440977575.0\n",
      "Iteration 11\n",
      "Iteration 12\n",
      "Iteration 13\n",
      "Iteration 14\n",
      "Iteration 15\n",
      "Iteration 16\n",
      "Iteration 17\n",
      "Iteration 18\n",
      "Iteration 19\n",
      "Iteration 20\n",
      "16.465544652633543 1.4942736290572998 [3.63564195e+08 1.41784559e+08 1.02854522e+08 1.13898111e+08\n",
      " 2.08399631e+07 1.38488876e+07 5.54145514e+06 7.39854145e+06\n",
      " 3.02391563e+06 1.21892687e+06 4.90150119e+05 2.22095466e+05\n",
      " 1.30197006e+05 6.08268087e+04 2.63639533e+04 1.27042786e+04\n",
      " 6.73337260e+03 3.31555813e+03 2.10636333e+03 1.36269182e+03\n",
      " 0.00000000e+00 9.99999970e+09] 3352826690825825.5\n",
      "Iteration 21\n",
      "Iteration 22\n",
      "Success!\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import minimize_scalar\n",
    "from functools import lru_cache\n",
    "\n",
    "method, xopt = \"bounded\", \"xatol\"\n",
    "\n",
    "new_hist = hist[:-1]\n",
    "coverage_range = np.arange(5, len(new_hist), dtype=np.int)\n",
    "P = new_hist[coverage_range] * coverage_range\n",
    "\n",
    "@lru_cache(maxsize=None)\n",
    "def nbinom_pmf_range(lambda_: int, rho: int, bin_id: int):\n",
    "    stacked = np.zeros(len(coverage_range), dtype=np.float64)\n",
    "    lambda_ /= 100  # 2-digit precision\n",
    "    rho /= 100      # 2-digit precision\n",
    "    n = lambda_ / (rho - 1)\n",
    "    p = 1 / rho\n",
    "    start, end = bins[bin_id]\n",
    "    for i in range(start, end + 1):\n",
    "        stacked += nbinom.pmf(coverage_range, n * i, p)\n",
    "    return stacked\n",
    "\n",
    "def generative_model(G, lambda_, rho):\n",
    "    stacked = np.zeros(len(coverage_range), dtype=np.float64)\n",
    "    lambda_ = int(round(lambda_ * 100))\n",
    "    rho = int(round(rho * 100))\n",
    "    for bin_id, g in enumerate(G):\n",
    "        stacked += g * nbinom_pmf_range(lambda_, rho, bin_id)\n",
    "    stacked *= coverage_range\n",
    "    return stacked\n",
    "\n",
    "def func(lambda_, rho, G):\n",
    "    stacked = generative_model(G, lambda_, rho)\n",
    "    return np.sum((P - stacked) ** 2)\n",
    "\n",
    "def optimize_func(lambda_, rho, G):\n",
    "    # Iterate over all G\n",
    "    for i, g in enumerate(G):\n",
    "        G_i = optimize_func_Gi(lambda_, rho, G, i)\n",
    "        if not 1 < G_i < 9.9e9:   # Optimizer did not optimize this G_i\n",
    "            break\n",
    "    # Also remove the last bin since it is subject to marginal effect\n",
    "    G[i - 1] = 0\n",
    "    lambda_ = optimize_func_lambda_(lambda_, rho, G)\n",
    "    rho = optimize_func_rho(lambda_, rho, G)\n",
    "    score = func(lambda_, rho, G)\n",
    "    return lambda_, rho, G, score\n",
    "\n",
    "def optimize_func_lambda_(lambda_, rho, G):\n",
    "    def f(arg):\n",
    "        return func(arg, rho, G)\n",
    "    res = minimize_scalar(f, bounds=(5, 100), method=method, options={xopt: 0.01})\n",
    "    return res.x\n",
    "\n",
    "def optimize_func_rho(lambda_, rho, G):\n",
    "    def f(arg):\n",
    "        return func(lambda_, arg, G)    \n",
    "    res = minimize_scalar(f, bounds=(1.001, 5), method=method, options={xopt: 0.01})\n",
    "    return res.x\n",
    "\n",
    "def optimize_func_Gi(lambda_, rho, G, i):\n",
    "    # Iterate a single G_i\n",
    "    def f(arg):\n",
    "        G[i] = arg\n",
    "        return func(lambda_, rho, G)\n",
    "    res = minimize_scalar(f, bounds=(0, 1e10), method=method, options={xopt: 100})\n",
    "    return res.x\n",
    "\n",
    "G0 = np.zeros(len(bins))\n",
    "l0 = 17\n",
    "r0 = 1.6\n",
    "print(l0, r0, G0)\n",
    "def run_optimization():\n",
    "    ll, rr, GG = l0, r0, G0\n",
    "    prev_score = np.inf\n",
    "    for i in range(10000):\n",
    "        print(\"Iteration\", i);\n",
    "        ll, rr, GG, score = optimize_func(ll, rr, GG)\n",
    "        if score / prev_score > 0.9999:\n",
    "            break\n",
    "        prev_score = score\n",
    "        if i % 10 == 0:\n",
    "            print(ll, rr, GG, score)\n",
    "    print(\"Success!\")\n",
    "    # Remove bogus values that are close to 1e10\n",
    "    final_GG = [g for g in GG if 1 < g < 9.9e9]\n",
    "    return ll, rr, final_GG\n",
    "        \n",
    "ll, rr, GG = run_optimization()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 500,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[363562395.98472995, 141791180.85236096, 102839547.96888512, 113924375.8267978, 20802719.768104754, 13890524.572856257, 5507316.05049201, 7415375.966161654, 3022713.950959648, 1219027.1216367115, 490129.2270706757, 222020.93582455214, 130134.20417864512, 60770.60368314885, 26365.8733778964, 12775.758506410042, 6664.702551278448, 3383.907247071877, 2038.4307060221818, 1297.0701514985467] 16.465773621960782 1.4947515297228837\n",
      "Genome size=2830468755\n",
      "Copy 1: 363562396 (12.8%)\n",
      "Copy 2: 283582362 (10.0%)\n",
      "Copy 3: 308518644 (10.9%)\n",
      "Copy 4: 455697503 (16.1%)\n",
      "Copy 5: 104013599 (3.7%)\n",
      "Copy 6: 83343147 (2.9%)\n",
      "Copy 7: 38551212 (1.4%)\n",
      "Copy 8: 59323008 (2.1%)\n",
      "Copy 9-11: 90681419 (3.2%)\n",
      "Copy 12-16: 85331899 (3.0%)\n",
      "Copy 17-23: 68618092 (2.4%)\n",
      "Copy 24-32: 55949276 (2.0%)\n",
      "Copy 33-45: 65978042 (2.3%)\n",
      "Copy 46-64: 63505281 (2.2%)\n",
      "Copy 65-91: 55526529 (2.0%)\n",
      "Copy 92-128: 51997337 (1.8%)\n",
      "Copy 129-181: 54750531 (1.9%)\n",
      "Copy 182-256: 55580677 (2.0%)\n",
      "Copy 257-362: 66874796 (2.4%)\n",
      "Copy 363-512: 85120229 (3.0%)\n",
      "Copy 513+: 333962777 (11.8%)\n",
      "Ploidy: 4\n",
      "Repeat content: 50.1%\n",
      "SNP rate: 1.15%\n"
     ]
    }
   ],
   "source": [
    "print(GG, ll, rr)\n",
    "genome_size = int(round(total / ll))\n",
    "inferred_genome_size = 0\n",
    "for i, g in enumerate(GG):\n",
    "    start, end = bins[i]\n",
    "    mid = (start + end) / 2\n",
    "    inferred_genome_size += g * mid * (end - start + 1)\n",
    "inferred_genome_size = int(round(inferred_genome_size))\n",
    "genome_size = max(genome_size, inferred_genome_size)\n",
    "print(\"Genome size={}\".format(genome_size))\n",
    "copy_series = []\n",
    "for i, g in enumerate(GG):\n",
    "    start, end = bins[i]\n",
    "    mid = (start + end) / 2\n",
    "    copy_num = start if start == end else \"{}-{}\".format(start, end)\n",
    "    g_copies = int(round(g * mid * (end - start + 1)))\n",
    "    copy_series.append((mid, copy_num, g_copies, g))\n",
    "    print(\"Copy {}: {} ({:.1f}%)\".format(copy_num, g_copies, \n",
    "                                         g_copies * 100 / genome_size))\n",
    "if genome_size > inferred_genome_size:\n",
    "    g_copies = genome_size - inferred_genome_size\n",
    "    copy_num = \"{}+\".format(end + 1)\n",
    "    copy_series.append((end + 1, copy_num, g_copies, g_copies / (end + 1)))\n",
    "    print(\"Copy {}: {} ({:.1f}%)\".format(copy_num, g_copies, g_copies * 100 / genome_size))\n",
    "    \n",
    "# Determine ploidy\n",
    "def determine_ploidy(copy_series, threshold=0.15):\n",
    "    counts_so_far = 1\n",
    "    ploidy_so_far = 0\n",
    "    for mid, copy_num, g_copies, g in copy_series:\n",
    "        if g_copies / counts_so_far < threshold:\n",
    "            break\n",
    "        counts_so_far += g_copies\n",
    "        ploidy_so_far = mid\n",
    "    return int(ploidy_so_far)\n",
    "\n",
    "ploidy = determine_ploidy(copy_series)\n",
    "print(\"Ploidy: {}\".format(ploidy))\n",
    "\n",
    "# Repeat content\n",
    "def calc_repeats(copy_series, ploidy, genome_size):\n",
    "    unique = 0\n",
    "    for mid, copy_num, g_copies, g in copy_series:\n",
    "        if mid <= ploidy:\n",
    "            unique += g_copies\n",
    "        else:\n",
    "            break\n",
    "    return 1 - unique / genome_size\n",
    "\n",
    "repeats = calc_repeats(copy_series, ploidy, genome_size)\n",
    "print(\"Repeats: {:.1f}%\".format(repeats * 100))\n",
    "\n",
    "# SNP rate\n",
    "def calc_snp_rate(copy_series, ploidy, genome_size, K):\n",
    "    # We can calculate the SNP rate s, assuming K-mer of length K:\n",
    "    # The probability of kmer containing no SNP is (1-s)^K. The proportion of Kmers under the 'het' peak (cov 2 to cov5), counting only distinct Kmers L=200+400+600+400=1600. The proportion lacking any SNP is thus 1-L/2G=1-1600/(3150*2)=0.75. Note the factor of 2 before G is to account for both alleles. Finally, This ratio equals to (1-s)^K.\n",
    "    # s = 1-(1-L/G)^(1/K)\n",
    "    # L: # of unique K-mers under 'het' peak\n",
    "    # G: genome size\n",
    "    # K: K-mer length\n",
    "    L = 0\n",
    "    for mid, copy_num, g_copies, g in copy_series:\n",
    "        if mid < ploidy:\n",
    "            L += g\n",
    "        else:\n",
    "            break\n",
    "    return 1 - (1 - L / genome_size) ** (1 / K)\n",
    "\n",
    "K = 21\n",
    "snp_rate = calc_snp_rate(copy_series, ploidy, genome_size, K)\n",
    "print(\"SNP rate: {:.2f}%\".format(snp_rate * 100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 477,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 600000000.0)"
      ]
     },
     "execution_count": 477,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEXCAYAAABPkyhHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd5zcdbX4/9eZvr1vssmmbQoJhBJC6HBDUcoFRSniRYr3KoogTa96r3Ct+KVcvVwb4tULolcpFsoPEYgU6RggQArp2Wyyvcxsmz7v3x+f2WST7GZnd6fveT4e89jJZz4zc+aTmTnzeZfzFmMMSimlpjZbpgNQSimVeZoMlFJKaTJQSimlyUAppRSaDJRSSqHJQCmlFJoMlFJKkQXJQETqROR2EXleRPpExIjIykk8nk1EPi8i74pIv4i0iMjjInJMEsNWSqm8kvFkABwCfBWoB95LwuPdAdwTf6ybgf8GjgBeFpHDkvD4SimVdyTTM5BFpARwGWO6ROQC4E/AacaYFybwWDagF3jKGHPxsO1LgfeBbxtjvpGcyJVSKn9k/MzAGNNnjOkaa79488+XRWSDiATjzT8/EpHiYbs5gEKgbb+7t8b/+pMUtlJK5RVHpgMYh18C/wT8L3A3sBC4DjhURM40lpCIvA5cJSKvAX8DKoFvAy3ArzITulJKZbecSAYicgpwFXCRMeYPw7b/HXgQOAv4S3zzFcBDwG+GPcQm4GRjTEtaAlZKqRyT8WaiBF0EdAMvikj10AXrl38UWDls315gLfAj4OPAFwAP8ISIVKY1aqWUyhE5cWaA1SRUCXSMcnsNgIg4gL8Cq4wxNw3dKCKrgHXAl4CvpzZUpZTKPbmSDGxYbf5XjHJ7c/zvqcBS4PrhNxpjNovIBuCklEWolFI5LFeSwVbgNOAlY0zwIPtNi/+1j3Cbk9x5vUoplVYJ9xmIyAoReVJEeuIze98VkatSGNtwvwdcwNdGiMstIqXxf26K/710v32Oxprc9k4qg1RKqVyV0C9lETkHeAx4AbgVCAOLgFnJCEJEbolfXRL/e7mInAx4jTE/NsY8LyK/AL4pIsux+gVi8RguAS7D6id4S0SeBf5FRMqBVUAd8EVgAGs2slJKqf2MOQNZRMqwfnE/aIy5ISVBiIwWRKMxZm58HwE+B1yNlTSCwHbgSeBuY0xnfL8C4MtYZwfz4vu9BNxijElGuQullMo7iSSDa7AmedUaY3zx8hH9JtN1LJRSSiVNIn0GZwIfAOeKSBPWOP7ueKXRkTpqlVJK5ZhE+gwWYPUN3A/cidUJex5WpVEPcOP+dxAR7xiPWQYYrMSilFIqMaVAzBiT9JGRiTQTbQUagK8ZY+4Ytv1h4AJgxlB7/bDbEkkGlJWVTSRmpZSaknw+H4AxxiS9ekQi2WWo0ufv9tv+f8DFwLHAn4ffYIwpP9gDioi3rKyszOsdK2copZQaUl5ejs/nS0mLSiLZZai42/5loYf+XZG8cJRSSmVCIsngrfjfmfttr4//Ha1ekFJKqRyRSDJ4JP73X4Y2xMf8fwZrItfrKYhLKaVUGo3ZZxCf1fsA8G8iUgu8Dfwj1hoCXzHG6IggpZTKcYkOT/ossBO4Mn7ZBnzeGHNvqgJTSimVPmMOLU3Jk+poIqWUGrf4aCLfWCM2JyJXVjpTSimVQpoMlFJKaTJQSimlyUAppRSaDJRSSqHJQCmlFJoMlFJKoclAKaUUmgyUUkqhyUAppRSaDJRSSqHJQCmlFJoMlFJKoclAKaUUmgyUUkqhyUAppRSaDJRSSqHJQCmlFJoMlFJKoclAKaUUmgyUUkqhyUAppRQJJAMRWSkiZpTL4nQEqZRSKrUc49j3buCt/bY1JzEWpZRSGTKeZPCiMebRlEWilFIqY8bVZyAiJSIyngSilFIqB4wnGfwa6AX8IvKMiByeopiUUkqlWSK/8kPA74GngE7gCODLwMsissIYs2n/O4iId4zHLBtvoEoppVJHjDHjv5PIkcBq4GFjzGUj3D5mMigrK8PrHWs3pZRSQ8rLy/H5fD5jTHmyH3tC8wyMMe8Cq4AzRrm9/GAXwDeJmMetMxTimk2beLa7O51Pq5RSOWMyk86agMpkBZJKdzY18bPmZs59/33+1NGR6XCUUirrTCYZNABZ/81qjOGPbU0ARIzhkvXrebO3N8NRKaVUdklkBnLNCNtOBk4Dnk5FUMn0weAgW0PW9TK8RIzhd+3tmQ1KKaWyTCKjiR4SkUHgVazRREuBq+PXv5m60JLjD23bAailjTNZxW+5jJe6m4AFmQ1MKaWySCLNRI8CNcCXgJ8AFwK/BVYYY3amMLak+GPbDgBOlrdZ7vED8N5gjGAslsGolFIqu4x5ZmCM+SHwwzTEknTtwSDvBIsBOL+yhKWFh0IThLHzTl8fx5fpdAellII8L2G9pmfznuvnzD6bhopjmcFuAF7xZn3ft1JKpU1eJ4ONfdYoomo6qS09mtLSY1nCBgBe7s76Fi6llEqbvE4GWwatGc71Ni8iNhyOUo5yWfPd3uwPZTI0pZTKKnmdDLb6gwDMdob3bDu+tBSA5qiH1mAwI3EppVS2yetk0Bi2+scbPHv7yZdXLtlzfcNgf9pjUkqpbJTXyWBX1DoLWFC4d9TQtPLjqaUNgHW+xozEpZRS2SZvk0F3wIsXKxkcUlK3Z3tBQQP18dU6N/briCKllII8TgbrfXuHlS4p2zvbWMTGXIfVPLR5cDDtcSmlVDbK22SwsXcXAAX4qSuYts9tDW4BYFtI0h6XUkplo7xNBpsHewCYZevBZtv3ZS4qKAKgMVJEbAKL+yilVL7J22SwLT6sdI7zwPkEi0tqAQjhpEmHlyqlVP4mg6FhpXPdB5ZfWlzagI0oAB9oJ7JSSuVvMmiJWk1B8woKD7itvPgQptMKwDpfU1rjUkqpbJSXycCYGN3xYaV1npIDbnc6K6kX64xgY39XWmNTSqlslJfJoC/YzgBW6eoZnqoR95nnsNY22BLQPgOllMrLZLB7oGXP9ZmFtSPuM99jB2B7yJ6WmJRSKpvlZzLw7+0UHjUZxEtUNEeLMDq8VCk1xeVlMmj2W3MMXIQoczhH3KehqAaAIC7aQlrOWik1teVlMmgNWuUmKmUAkZFnGS8qnb3n+paB9rTEpZRS2Sovk0FbKABAtW30zuGaogZKsRa62dy3Oy1xKaVUtsrLZNAejgBQ7YiOuo/DUUEdVt/CloHOtMSllFLZKi+TQUfEahqqdY5eiE5EqLdbzUk7/Fq9VCk1teVlMuiMugCY5nQddL9ZrhgAjcHRzyCUUmoqmFAyEJGviIgRkTXJDigZumJWCYrp7gNLUQw3122NNGqKHDxpKKVUvht3MhCR6cAtwEDyw5m8SKSfbqw5BHWesoPuOy8+16AlVqqlrJVSU9pEzgxuB1bHL1lnINBCbzwZzCwYuRTFkAVF1oS0ME6ag/6Ux6aUUtlqXMlARI4FPgXcnJpwJq/Zv7cUxYzCaQfZExaVzt1zfVPfrlSFpJRSWS/hZCDW7K0fAb8yxhy0r0BEvAe7AAdvv5mEFv/eYaJ1nuKD7ltROJcKugHY0tecqpCUUirrjefM4ArgUKz+gqzV7PcCYCdKhePAhW2Gs9s91ImVDLbGl8lUSqmp6ODflnEiUoLVV3C7MaZlrP2NMeVjPF7Kzg5ag9acgUoZwDZKKYrh6h1+1odhR0DnGiilpq5EzwxuAULAD1IYS1LsmX1sT2ydgtlOaxRRoy5roJSawsY8MxCROuBG4FZg2rDCbx7AJSJzAZ8xJivaWboj1gSySntiE8nmetwwCLsinlSGpZRSWS2RZqJpgAu4I37Z3/b49q8lMa4J64nngIoE16xpKKyAbmg1ZUSNwZ5A01K+2h0Msqqnh+d7eogCJ5SWclFNDbUunZSnVL5LJBlsBz42wvbvAkXATcCmZAY1Gd6olQXG6jwesqB4OhAjioOdg17mFVWkMLrs9KrPxzWbNvHewL7zCH/T1sat27fz00WL+ETtyIsEKaXyw5jfmMYYH/Do/ttF5EYgYow54LZM8sasEhOVzpEXtdnforJ5wFYANvU1Trlk8FRXFxeuW4c/ZtVpKrHbWVlejlNgVY+X7kiES9ev562+Pu5oaBh1fQilVG7Lq0J1xhh6jRuAKufB6xINKfXMpIouADb3taYstmz0u7Y2PrJ2Lf5YjAaPh8fmDfC3ip/y1d4z+GLnEfxv9AJOsL0PwF1NTVy/ZYsuEapUnppwMjDGrDTGHJXMYCYrGu2nD2uiWbW7KKH7iNiYYbPmJmz3+1IWW7b5n+ZmLtuwgYgxHOox/Fi+ROn28/B2Pkg4bK3zUEU334ndyNk8BcCPd+/mnmadnKdUPkqsYT1HhMNd9FIKQLW7NOH7zXIEeT8E2wNTY3zpk11dfH7TJgywzNnOtwL/QgHW2g4VFR+muvpjFBUtxZgwPt/LfH33PQTCHl7gNG7c/AHHlRSxvPSgU0mUUjkmr5LBYKgDP1bzUI078bb/2S6BEDSF8qrVbETv9/dz6fr1xIAltkZuC3+eAgKUlp7IwoU/oqTk6H32r6g4jfr6G7lr079yQfsudlPPpe8+z/qTPorTlv/HS6mpIq8+zZ3B7j3Xq92JT3CeV1AAwO5oQdJjyiYD0SgXr1tHfzRKLR18J3YzhRJh/vzvs2zZSwckgiEORwnLl9zDj6c1IcTYEq3gjvW/TXP0SqlUyqtk0BHY2+ZfleBoIoCGwkoA2kw5oWj+rnp23ebNbPT7cRDm29zCdKedI498nlmzbkbk4G8FEeGji7/Gxz3WKOK7OkvZ0fVKOsJWSqVBXiWDzmDvnuuJzjMAWFgyE4AYdrYP5OeIogfb2ri/1XptV/NzjnAHWLbsZcrLT074MUSEHx75SYrw00spt274M7FYKFUhK6XSKL+SQcgqNuchhMee4BRkYGFpA4I1zn5TX2NKYsuknYEAn9u4DoDjeJ0rCtaxbNkrFBYuGvdjzSgo4+YZVufx7yMn817jj5Iaq1IqM/IqGXSHAwCU2sb3a7XQVUHN0LoG/R1JjyuTosZw6fuv0BuzUU4P3yr8C0cf/RIeT/2EH/NLDSsolRABCvivxk0EgzrcVKlcl2fJIAxAuS087vvOtFlNTNv9vWPsmVu+velFXhuwzpJudT/E6cv+gNN58OVAx1LmcHB9/WwA/sD5rN1+56TjVEplVl4lg554xdIK+/hnydY7rQTSGIwkNaZM+lvrm9zWYr2eC+3Pc+3R/zXpRDDkptnzKZQoAxRzf2sLgcDOpDyuUioz8ioZeGNW3ZwKx/jr58xxW4eiKZR4X0M28w02ctXGjURxMEd28/Ojr8Ltrkva41c6nVw+zXq8P3E+23d8N2mPrZRKv/xKBvGKpeXjGEk0pKHAKmOxK3rwdZNzQSTSy43v/IztZhZ2ovx6yeFUFs1L+vN8cdYcAHZTzxOt6wiF2pL+HEqp9MirZOCLVyytcrrHfd+GomoAOinHH8nd4ZKxWIRH3ruBB8JnAnBDrXBK7REpea7Dioo4rcwq+/FHzmf37ntS8jxKqdTLm2QQjQ7Si1Wcrso5/pnEi0pmAWCwsaVvRzJDS6sNW77Ev/eeSgw7i10hvrd4ZUqf74v11nH7OytYs/tBotFASp9PKZUaeZMMwuEu+igBoNo9/qae+aXzsGF1QG/ua0pqbOnS2voA32seZAfzsBPj10uPx53i+kHnVVVR47ATw86TkWNob/9dSp9PKZUaeZUMhiqWVo2jYukQl93NtPhcg639XUmNLR36+tbw0MYf8zs+CcBXZs/hmNLxH4fxctpsfGq61ZH8F85m9+57U/6cSqnky5tkEAr30B9fy6BmHEXqhptpt8o4bw/0Jy2udAiHe3hl7ZV81/wrBhvLigr4xtzkdxiP5tPTpwNWR/Jr/QP097+ftudWSiVH3iQDb6iHGNZoool0IAPMclnNRI3B3ClWZ0yMDRuu4DvBS+ikhiIbPHTY4SlvHhru8OJilhdbifgpzqGl5X/S9txKqeTIm2TQMaxIXeUEhpYCzHVbo5Gawq6kxJQOjY3f45fdDl7lJAB+umgxCwsTW/Izmf65zmoqep7T2N76CNGoP+0xKKUmLm+SQVdob9NO5TjKVw+3sMgqwLYjVkksvkB8NuvufppVO+7nHq4B4FPTpnFFvMkm3T5ZW4tbhAAFPBc9kq6uJzMSh1JqYvIoGVhDGu3EKBlHxdLhjiqfC0AfJezs35GkyFLD79/B6nX/zLf5D8K4mO9x85OFCzMWT4XTyQXV1lyNpzhHRxUplWPyJhnsqVgqQUTGX44C4IjyRXuGl77TsylpsSVbNBpg3boL+X70CpqYjVPgocOWUjrB5rFk+XS8qeg9juT9zreIRHxj3EMplS3yJxnEi9SVTaBi6ZAiZwEzxBpW+l5vdpZWMMawefO1/KG/mmc4C4C75i9geUlJhiODMysqqHNZTXTPcRIdHX/KcERKqUSNmQxE5BgR+ZOINIqIX0RaReQvInJiOgJMVE/EauOvsE+urX+hcwCA9YPZ2QHa0vIL3mx9lru5EYCPVFVx/cyZGY7KYhfhkppaAJ7jdNrbdZ1kpXJFImcG8wEH8D/AdcBdQC3wNxH5UApjG5eeqNU0VDbJoqNLCqwH2BzMvhFFXu/LrN70b3yd2whQwCy3m/sWL55ws1gqfKLWSgbbmM+ans1avE6pHDFmI7Mx5iHgoeHbROQeYBtwA/BsakIbH2/M+hKvdEwuGxxWXAU+2BqrJhaLYLNNrB2+J9jHrxtf5PGubgLGTqXDxhnVc7lk5tHUucc/D2JwcCNvvX8Rt3Iru5iFxyY8cthhEx45lSrHl5Yyx+2mMRjkef6BM9ofpr7+i5kOSyk1hgn1GRhjBoEOoDy54Uxcb9T6Upzsl+PQiKJeSmns2zqhx3il/S0WvbaKG5qL+WtwNq+EZvLEYB037gwy+7WXuGXjS4THMXTV79/Om2vO49+jN/EeRwLwq8VLOC4N5SbGS0T2nB08x+m0temoIqVyQcLJQERKRKRaRA4Rke8BS4G/pi60xBkTxYdVqbTS6ZnUYx1ZsRDB+qJ+t2fzuO//p8an+dD6NjqpwE2As5zr+ZfC9zjT/hbF9BHBwW0tUZa++ihPt4+dbAYHN/HkOxdxQ+g63uQ4AH4wfz6XxL9ws9FQMmhiNu/0teP3b89wREqpsYynDeQ+4ML49RDwM+B7I+0oIt4xHmtixYNGEYl491YsdY2/fPVwRQ4XM6Sb3aaad3tbuWAc9327ay2Xb4/gp4hpdPHEYfNYUXM2YI0Cau5+ka988CK/C5/CpkglZ69v4sTt6/haw3H8Y3UttmFt/8YYmtoe4baNf+EB810C8WT3/fnzuWnWrEm9xlRbVlzMwoICNvv9PM9pfKj9IebM+Vqmw1JKHcR4ksG3gHuBeuBywA04gWAK4hqXcLh7TzKock1+iOWhLj+7g7C6P/ERRZ0BHx9du4EBaqiih5ePPpYFpXu/tEWEmVUr+c2Jp/Cp7b/hKzv9rGUxr/pL+ci6DTQ43+YfK4qodroJhVpY793CS+E5dHIFAHVOuPeQpZwfn9iVzYaair7b2MhznE5Hx39pMlAqyyWcDIwx7wPvA4jIb4DVwP3ARSPse9C+hPiZQ9LODiKRnmFrGUy+Hf3EshKebYe3QjUYE0Vk7E7pm9b+iV1mLi6CPLKoZp9EMJyInXMaruSM+i4e+OAn3NtdyGqOYVu4gB+1xwA/VlfMMQA4iHJtXQ3fmn8oZRmeVDYel8aTQSt1/L2/j6WBJjye7D6jUWoqm2gHchh4DPi4iEyuXSYJeoPdBLH6Cqpdk1/D+PSaJQC0MJ1N3WvG3P+N9jX8tt/6ovvXit2cNuPkMe/jclXxmSP+g+ePOZsnav/Kx23PsII3OYQPWMImTnZu45t1hg+OPYG7DzkipxIBWEtiLo0XzHuO0+nsfDTDESmlDmYy3zAFgAAlWD9nM6Yr6IN4MqhyTX5+wPGVC3CzgyBunmtfxyFVy0fdNxaLcd3Gt4nRwBxp5ZbDPjmu5youXsp5hy7lPCAa9ROJeHG5piGS+5PDP1Fby9odO3iBlbR3/FSHmCqVxRKZgVwzwrZS4GKgyRjTnorAxqMz1LfnekUSfkG77HaOdHYC8LLv4H3h929/ktXRBgDumluFxzHxEyW7vQC3uy4vEgHsHVXUSQ0v+3oIhTozHJFSajSJfOs8JCJ/FpFbROQzIvItYC0wG/hyasNLTGdocM/1ZCQDgBNKrDOMvwdLMcaMuM9A2M/Xd1n956c6t3LxnKyZkJ0VFhYWsry4CIDnWElX1+MZjkgpNZpEksFvgELgeuAe4AvAu8BpxpiHUxhbwnriFUuLCOJM0gpfp1dZv/a3mNm09K4bcZ9vbniUVlONgzA/XHxsUp4333yidhoAL/IPtHZov4FS2WrMb05jzP8aY1YaY2qNMU5jTI0x5nxjzIvpCDAR3WGrUmmZLZS0x1xZezguwhhs3LfjuQNu39a3mx93VwBwZdEWjqw6PGnPnU+GJsd5qeC57nYikb4x7qGUyoS8aJzujkQAKLcnb+3iUqeL84p7AHigp5hoNLDP7deve4YAHirp4fbDPp605803czweji8Zaio6he7upzIckVJqJHmRDHriOaDcPnLb/kRdP9caRbSJBp7dtXcZx0can+fJwDwAbpk2SHXhtKQ+b7755DRr0ZuXOIXmdm0qUiob5UUy8MXLV1fYk1vK+dSqBg6xdwBw186dBEM9vNG5nqu2WyNpD7ft4PpDPpHU58xHF9fUIBj6KOWZ7pYDzrKUUpmXF8nAG7MqlVY4kzsxS0T4wox4Bc7oMk547WHOXLudQQqppZM/Hnky9gmWuJ5K6txuTi21mopWxU7A682K+oZKqWHyIhn4Ytb6AFXO5C9Ic23DqVxRbq1+9o45hH6KKKGPPy6ezoKyBUl/vnz1yen1ALzCSTS1P5bhaJRS+8v5ZBCN+unF+tVZ6Ux+ZQy7CPcfeS7/MT1Mg72bGytaWb/8CE6aPnbJCbXXhdXV2DEMUsSfO3cTi0UyHZJSapicb+MYXqSuylWUkucQEb61+EN8a3FKHn5KqHa5OL2skGd9fp6JruAa38tUVKzMdFhKqbicPzMYXr66JgkVS1Xq/NP02QC8zvE0alORUlkl55NBMNxN/57y1VmzCqcawQXV1biIEcTDYx27Ri3zoZRKv5xPBl3BvYXkql2TW/JSpVa508mHyq1+nWcjR9LXtzrDESmlhuR8MugI9u65nqwidSp1Lquzaj69ybFsa9fCdUpli5xPBt2h/j3XK53ODEaiEnF+VRUeiRLGxR/admpTkVJZIueTQVfYmg3sJEJhkiqWqtQpdjg4p9yaF/JM+FAGB9dnOCKlFORDMghZ6wmUSRCR5JajUKlxWd0iAFZzDBtbn8hwNEopyINk0B2xqtSV2XUSU644t6qKIokQw87v23dmOhylFHmQDHoiMQDKbckrX61Sq8Bu5x/Lrf6dvwQX4Pdvz3BESqmcTwbeeA6o0IFEOeXymYcCsIajWKdNRUplXM4nA1/MDkCFw57hSNR4fLiyijIJYrDxcFtjpsNRasrL+WTgjVmVSisdOqw0l7hsNs4vt07n/hyYRyCwK8MRKTW15XQyMCaGz1izjiud7gxHo8brivojAFjHUl7brbWKlMqknE4G0Wgf/RQDqatYqlLnjMpq6u3WpMH7WjszHI1SU1tOJ4NwuJterEql1a7iDEejxssmwpU11v/bE+Ej8A1szXBESk1dYyYDEVkhIj8RkfUiMiAiO0XkQRHJ+DJf+5Sv9pRlOBo1EZ+feyw2Ynip4HeNqzIdjlJTViJnBl8FPg6sAm4Afg6sBN4RkSWpC21s3mAPEayO42qXrmWQi+o9hfyDpw2AB7p0rohSmZJIMvgBMMcYc70x5hfGmO8CpwBOrESRMZ0h357rWqQud32mbhYAr0cX80H32gxHo9TUNGYyMMa8aowJ7bdtM7AOyOiZQVeob891TQa566L6E6jCi8HGzxr/nulwlJqSJtSBLFZFuGlARoeAdAYHrXiIUaZrGeQsl93ORSVdADzYW0YkFstwREpNPRMdTXQZMBN4eKQbRcR7sAuQlN7e7nj56hIJYteKpTntc7OtOQdtppLHm3UFNKXSbdzJQEQWAz8BXgZ+nfSIxqErbLVelUk4k2GoJDiq+hiOtm0E4Oe7Nmc4GqWmnnElAxGZDjwJ9AAXG2NGPJ83xpQf7AL4RrrfePVErLLVZXYdhZLrRIQrq62mvlWBabQG+se4h1IqmRJOBiJSBjyF1cRzljGmNWVRJagnYi2ZWKE16vLClfPOpoReojj42bYXMh2OUlNKQslARDzAE8Ai4DxjzMaURpWg7qiVBaqd2l+QD8oKZnK+x5qFfF9niJiuj6xU2iQyA9kOPAScgNU09HrKo0pQd8wqTletw0rzxudmzgdgZ6ySR1uz4jeHUlNCImcG3wc+gtVEVCkinxp2uSC14Y3OmBheUwBAjasgU2GoJDt55jkcK2sAuLNRk4FS6ZLI4Pyj4n/Pj1+GawQeTWpECYpEvPjiRepqtWJp3rDZnFxTFeDNTngjUMbfe32sKNW6U0qlWiIzkFcaY2SUy9w0xDiicLgLL+UA1Lq1LlE+uWjueTRg9R3cvu2dDEej1NSQsyWsB0Md9Mcrlk7zVGQ4GpVMxcVLucLzFgCPeWM0BgIZjkip/JezyaA90LPn+nRPSQYjUalw1ewVVNNBFBs/0L4DpVIud5NBcO+8NR1NlH/qp13KRbanAfhlawfesM4yVyqVcjgZWBVLhRiVWqQu79jtBXxuRh0FDDJgHNy9a2emQ1Iqr+VuMghZ7cilEsBhy9mXoQ5icf1n+Vh8sNoPmhrp1rMDpVImZ79FO0NWkboKW2iMPVWu8nhmc21VkCL66YvZ+M+mpkyHpFTeyt1kELGK01Xa9ddiPjt87o1cxO8B+OGuRjpCmvyVSoWcTQZdEaseUZUWqctrJSXL+Ex5FyX0MhAT7tSzA6VSImeTwVCRuipnzr4ElaDD5tzIJfF1lH6yayetwWCGI1Iq/+TsN+ReOtsAABXySURBVGlPzAVAtUOHlea78vLT+XRpC2V48Rvh242NmQ5JqbyTk8nAGEOPKQS0SN1UICIcOu8WLuP/ALi3eTfv9eviN0olU04mg2i0f2+ROrcWqZsKKipO46qybmbTSAzh+s2bMbregVJJk5PJIBTqwIdVyXKaRytaThWLGm7jOn4MwIs+Hw93dGQ4IqXyR04mA2+wkxDWwjbTPJUZjkalS1nZ8ZxbM5+TeQmAL2/ZwkBU179WKhlyMhm0+vf+Ipzu1jODqaSh4f/xBX6BkxC7QiH+n3YmK5UUOZkMWvzde67XuFwZjESlW0HBPI6fcwWX8iAAdzbt5H3tTFZq0nIyGewOWh/+QoKUaJG6KWf27H/nn92vUU8TYQNXfLCBUCyW6bCUymk5mQxaglaRuhrbQIYjUZlgtxdwxCF381XuwEaUNf0DfHPHjkyHpVROy81kELY6DWsdkQxHojKlsvIszpx+PJ/gIQBu37mTp7u7x7iXUmo0OZkMWsNWXaLpDh1nPpUtWPADrnH9lUNZhwEu37CeZi1VodSE5GQyaI9ancYzXNpfMJU5HKUcfuj93MptlNBLRzjCP61fT0T7D5Qat5xMBh0xa9ZxnVtLUUx15eWncGLDjXyVOwBrMpr2Hyg1fgklAxGpE5HbReR5EekTESMiK1Mc24iiUT+dVABQ7ynJRAgqy8ya9WU+Wl3LRTwCwG07d/J/bW0Zjkqp3JLomcEhwFeBeuC91IUztoFAK954MphVqLOPFYjYWLz419xc9AZH8Q4An/5gAy/09GQ4MqVyR6LJ4C2g2hizELgrhfGMaZe/dc/1+sJpGYxEZROHo5hlhz/K7a57mcMOwgY+tvY91g/o8GOlEpFQMjDG9BljulIdTCKaBveGMaugPIORqGzj8dRzylGP833HXVTShTdqOOfd1ewKBDIdmlJZL+c6kHf5ewEoxK+zj9UBCgsXcdZRv+NOx5148LMzZDj17ddo0oSg1EGlJBmIiPdgF2DC1eWag9Zpf41N69GokRUXH8Gly37FHY67cRNge0g4efWLNPr9mQ5NqayVc2cGzaEwALV2nVykRldUtJirj/kld3vus84QIm5O+vsqtg54Mx2aUlkpJcnAGFN+sAvgm+hjt1q5gGkOnVikDs7jqeezK+7nZ+Wr8OBnd6yE41e/yCttqzMdmlJZJ+fODNojVj/BdKdkOBKVC+z2Aq448j/5zYxtFNFPpynjwxva+Om6e4jFtLaVUkNyLxnErFnHM9zuDEeicoWIcOGia1m1pIpa6WaQIq7tWMIlr95BW8/LmQ5PqayQU8kgFovSbKyJZnN0WKkap+OnncA7x57Bie52AP4QOYmT3v2Ax9+7jmBwd4ajUyqzEk4GInKLiNwCXBzfdHl823WpCe1AuwZ2MYhVl2hJyfR0Pa3KIzMKynjp+Iu5tc6OnShbWcBF3R/hmte/wfot/0Y4nBXTaZRKOzEmsTLQIjLajo3GmLnjelIRb1lZWZnXO76RHU/vfomzN1trGXhPXE6ZS2sTqYl7ydvD5etW0xh2AjCHHXzFdg/nzDqL+vobcTq13InKLuXl5fh8Pl98IE5SJXxmYIyRUS5zkx3UaD7o7wSgih5NBGrSTimvYP3xp/PlmdOwY2hkLtfG7uDzjWEef/UoNm++Hr9/W6bDVCotcqrPYPOgNeFstl0LkKnkKLTbuWvhEv6+/BiWF1uDE57hLD5pfsHNu5386Y1TWbfuYnp738hwpEqlVk4lg21Bq4lorlNLC6jkWlZSwhvLj+XnixZR63QQxsX/x/lczq+4vmMx9719NatXH0tz8/8QifRlOlylki6nksH2sDWctMGdU2GrHGEX4bMzZrDj+BP44YIFzHK7iGHnOc7gRv6bi/q/wG2bnuHRV49k48ar6e1dTaJ9bkplu4Q7kJP6pBPoQDbGUPzi0wzi4WfTd/K5xVekMEKlIBSL8du2Nu5pbubNvn3PBuazhRN5ldPdLZw6/STqpn+SgoL5GYpUTRWp7EDOmWTQEgwy47XXAHh+/iArZ52bqvCUOsCavj7ubWnh4fY2uiPRfW6rpIujWMNxbi8frl3CcTPOpaBgbmYCVXlNkwHwXMc2zli3E4CWZZVMLzsiVeEpNapILMZrvb081tnJ4x272Rw88PNTQi9LbM0cVeTkhKr5nFKznIbCYkS0hIqaHE0GwE+3v8K1jWGq6KT15A/jcJSmMEKlErN5cJC/dHfxbMdWXu4L0BPzjLhfmfg53BNiRWkFx1cuYEVpBXM8HmyaINQ4pDIZ5MzqMBv6vUAR9bRpIlBZY2FhIQsLC/li/SxixrBhoI+X2t/l9e7tvDMYZVNsOgEK8JkCXvYX8LI/Bm2bAHARZZYzyrwCD/MLq5hfWMI8j4d5BQXM83iocDj0bEKlTc4kgzcHrHbaw5wdGY5EqZHZRDisuJTDik/h8w2nADAY2MWbrc/zek8j7/QHWR+tZQsLCFBACDtbw3a2hmOs6u0A9n1vl9hizHHbmecpoKGwlPkFxXuSxVyPhyK7PQOvUuWrnEgGwViMNQGrJtGxnsEMR6NU4go99aycezkr51r/DoXa6Pa+yrvdb7Kpr51tgQF2RYtoZTot1NHKdAYoBqAvZmOt37DWPwg9B77vq2wBZjmCzHYZ5rodzPMUML+ojPlF1TQUTcfjGLnJSqmR5EQyWN3XRwjrV9BJ5VovRuUul2sa02s/xvTaj3FWfFs43MXAwDoGBtbi979Iu7+Fbf5+dgSj7IoW7UkSQ3/DuADoinnoCnlYEwL2rALbD/RjYyvVdDPT1s1Mez+zHCFmu2LMczuZV1DGnOI6Cj2z8Xhm43BUaHOUyo1k8JLXKj8xi53MLTs0w9EolVxOZxXl5adSXn4qAAuAE+O3RSL9hMNthELWJRjazm5/D9v8fhpDEXaG7DRFCtgdK6HZVNNBDTHsxLDTTg3tsRreiQFhYNgS0HYi1LCBabxAHZ3McgaY7YR5BW7mFZYzp6ieksIGPJ4GnM5qTRZTQE4kg791NwOwlLUUF1+b4WiUSh+HoxiHo3ifCW31wHEj7GtMlMFgN9sHWtk62M02fx/b/QEaQzEraUQL6IoVAhDFQSt1tFLHu2AlizAwCHSBjSiV7KSat5lOFzMdIWa5bcz1FNJQVMn84plMK5qHxzMXu12bo/JB1ieDmDG81hcAbBxl24LbPTvTISmVlUTsFHlqWOqpYWnVyPsMRqPsDARoDATYNtjFtoFOtvn7aQyG2Rmy0xEfGhvDTic1dFLDBwCR+GUAiC/5UMhWanmNOpuPGY4wdS4Hde5i6gsqqS+czuySOdQX1uHRju6ckPXJ4M3eXrwxqxbRcYVGT1eVmoRCu53FRUUsLiqCqipg0T63+6NRdgaD7A4GaRzsYdtAO41+L42BIE1hoTlSQBBr/YdBitjBPHbEgFD8sqfvIgxsAbZQygDV9gC19hi1Lid17kJmeMqZWVDBjIIKapwualwuqp1OHSGVQVmfDL7T2AhYtWCWlM7KcDRK5bcCu51DCgs5pLAQKiqAhn1uN8bQHg6zMzDIlr7meLLwsTMQpC0SoyPipNMUE2TvGuW9FNEbLWJblGEJYzB+2Xe5UQ8RKuwhquwRKu2GKoeNaqeDaqeLKpeHalcx1e5iql0lVLuKqXI6KXU4dPJeEmR1MnjN5+PP3d0AXMX9lJR8OsMRKTW1iQjTXC6muVysKC0HDhzQYYzBG/TR2L+dpv5d7PZ30Bzw0Rr00xqO0RF10k053VTipZxI/EwDIICDlqiDlugBD4t1ttETv+xlI0YJg5SJnzJbgHJbmHJbhHJ7jAoHVNhtVDjsVDqdVDvdVDoLqHYXUeEqwu0oxW4vjl9KsNmcIz3xlJCVyeDr27bx85YWusJhABaxkZN4haKiuzMcmVJqLCJChaecCs8yjqpedsDtxsQIhdoIh9sJBtvpCTbTEuihPdBLe2iQrnCUzij0RO10x1z0xDz4jIdeSumjhF5KiQ776ophw0cxPlMMUaxLQgIU0E0x/cMug5RIkBIJU2KPUWSLUWITSuxCid1GqcNJqd1JqcNNubuYClcFZe5q3K4anM5qnM7qnE0oWZkMwsbQGU8EguFqfo5NXBQVHZbhyJRSkyViw+2uw+2uo7gYqrCG0x6MMYZodIBo1Ec47MUX7qcr1E9naJCuUIDucIjuSJiuSAxvxNATBW/UjjfmwBdz4TNufKaQMPt+UfspxE8hHdQOe7L4JZbg6yFGATsp5AMKGaSIICW2MCW2GGV2KLXbKXa4KHG4KXUWUuosptRZQpmrjDJ3BZWuUma63ZRnuPxIViaDS2trObq4mBqXi9iuf8PZ9TYVFWdjtxdkOjSlVAaIyJ5htm73TIqBmeN8DGMM/liMrnCYrnCYnvAA3aEBekID9IT99ISD9IRDeKMRfJEo/dEYvTEYiAp9MTsDxs6AcWLY9wvbYGOQIgYp2rsxFr9ERosmhFV+ZN8SJA6iOInhkhguG7hEcIkNt81Ggd3BQDTh055xy8pkcHRJCUeXlBCLhXl13UNEgKqqj2Y6LKVUDhMRCu12Cu12Znk8QMm4HyNmDIPRKH3RKL3RKH2RCL5ImJ6QD2/Qizfchy80gDfspycSwhuJ0h+NMhCDgZgwaOwMGhd+PAxSSIx9R09FsBPBjt8wQnNXDFJYZTork8EQn+9vRCJWmevq6o9kOBql1FRnE6HY4aDY4aBun1tGmdgxgqEmr1Cog4FQJ+2BLnYFfPhCffgjffgj/fgjAwQiAwQigwSiAQLRIAGc3EeYcLJfVFxWJ4POzscAKClZgds9I8PRKKXU5A1v8iosnEcNMFZvqJVAennI8SQ+QimJK6GV5UXELSJ3iEiziPhF5HUROSMlEcVFIj7a2x8CoLpam4iUUlOXlUDKSPAre0ISfeT7gZuA3wA3YHWNPCUiJ6QoLnbs+CbhcDs2WyHTpl2ZqqdRSilFAs1EInIscClwkzHm7vi2B4C1wB3AqckOqr//fXbt+hEAc+bcisdTn+ynUEopNUwiZwYXYU39+8XQBmNMAPglcLKI1I12x4nasuUmIEpBwSJmzbo52Q+vlFJqP4l0IC8DPjDG9O+3/U1AgKOAluE3iMhYK92X+Xw+ystHW9M5RjRqx2bbjUjtKPsopdTU4vP5AFKyCHwiyaCO/atJWYYSwISH+fjir2x0AxN96KmmLP53jOOpEqTHM3n0WCZXGZCSacqJJIMCIDjC9sCw2/dhjBntJz+w98xhrP1UYvR4Jpcez+TRY5lcCbS6TFgifQZ+GFaPdi/PsNuVUkrlsESSQQswUifx0Lbm5IWjlFIqExJJBmuAxSJSvN/2oWVY301uSEoppdItkWTwe8AJfGZog4i4gU8Drxhj9MxAKaVy3JgdyMaYN0TkEeDO+JyCrcCVwBzgqtSGp5RSKh0SLVR3BfCd+N8K4D3gXGPMK6kKTCmlVPqISWF97FGfVIebJZUez+TS45k8eiyTK5XHMyPJQCmlVHZJXT1UpZRSOUOTgVJKKU0GSimlNBkopZQizckgE8tn5gMRWSkiZpTL4v32PVFEXhaRQRFpFZH/FpHCTMWeaSJSJyK3i8jzItIXP2YrR9n3IyLytogERGSniHxDRA4Yfi0i5SLycxHpEJEBEXlORI5K+YvJsESPpYjsGOW9evsI+07VY7lCRH4iIuvjr3uniDwoIgtG2Dehz/Rkv18TnWeQLPcDFwJ3A1uwJq09JSL/YIx5Lc2x5KK7gbf227ZnBnj8Q/RXYB1wM1APfBloAM5PU4zZ5hDgq1jvt/eAE0faSUTOAR4FngO+CBwO/AdQHf/30H424Mn47f8JdAFfAF4QkeXGmK0peyWZl9CxjHsL6/063Nrh/5jix/KrwEnAI1jHcjpwHfCOiBxrjNkA4/5M389kvl+NMWm5AMcCBrhx2DZPPOi/pSuOXLwAK+PH7oIx9vszsAsoHrbtM/H7np7p15GhY1cCVMWvXxA/FitH2G8d1heYfdi27wJRYOGwbZfs/38B1AA9wAOZfr1Zcix3AI8m8HhT+VieCLj227YQa2mA+4dtS+gznYzv13Q2E6V9+cx8JCIlozRdlAIfwvoQDV+V7gGgH+uDN+UYY/qMMV0H20dEDgUOBe41xkSH3fRTrKbUC4dtuwjrbOyxYc/RATwMXCAizmTFnm0SOZbDxZstDtZEOZWP5avGmNB+2zZj/ShZAuP+TE/6+zWdySCR5TPVwf0a6AX8IvKMiBw+7LbDsZr9Vg+/Q/wNtwbr+KuRDR2b/Y9dM9avsmX77fuWif/0GuZNrF/OB7T5TlEfxlqqcEBEtorI1SPso8dyGBERYBrQGd80ns/0pL9f05kM6thvreS4SS+fOQWEsKrH3gB8FPgW1mnhyyKyKL7PUOYf7Rjr8R3deI6dvo/H9h7wDawzqs9ifbndKyJf228/PZb7ugyYiXVmBGl+X6azA3ncy2cqizHmVeDVYZseF5EnsH4xfAPrTTR0/EY7xnp8RzfWsSvcb199Hx+EMeYjw/8tIvcBLwO3isg9xpih9ZD1WMbFRwX+BOs4/Tq+eTyf6Ukfy3SeGejymUlkjHkXWAUMDR0bOn6jHWM9vqMbz7HT9/E4xfth7sZKqicMu0mPJSAi07FGVfUAFxtjYvGb0vq+TGcy0OUzk68JqIxfHzodHO0Y6/Ed3XiOnb6PJ6Yp/rdy2LYpfyxFpAx4CigDzjLGtA67Oa3vy3QmA10+M/kagI749bVABDhm+A4i4sLqPFqT3tByytCx2f/YzcAa171mv32Xxzv7hjsOa4THllQFmeMa4n87hm2b0sdSRDzAE8Ai4DxjzMb9dhnPZ3rS36/pTAa6fOYEiUjNCNtOBk4DngaIt8OuAi7f7w1xOVCMNblFjcAYsw74ALhaROzDbroGiAF/GLbt91idcR8d2iAi1cDFwGPGmHDqI85eIlIZn0w2fJsH+FegDxg++WnKHsv4++whrGazi40xr++/zzg/05P+fk3regYi8jDWZJX/Yu/ymSuA04yumjYqEXkOGMTqRO4ElgJXAz5ghTFmZ3y/o+P7rMUab1wPfAl43hhzbgZCzwoickv86hLgn4D/BbYDXmPMj+P7nAc8jjUD+SGsY3wd1tyDLwx7LDtWJ99hWLNmO7Fmzc4Clhtj8v3X7EGPpYhcBXwd68tpB1CF9TlfBFxjjPnZsMeassdSRO7GGh34BHtHDw3pN8Y8Gt8v4c/0pL9f0zzrzgPchdW+FcAaA3tmOmPIxQtwPfAG1nT9MLAb60M4e4R9TwZeweowagN+CBRl+jVk+PiZUS479tvvAuCd+HuzCWsIr2OEx6uIfzA7scbSPw8cnenXmQ3HElge/4LbhTW6pRd4AasZZKTHm5LHMn5MEn1fJvSZnuz3q650ppRSSktYK6WU0mSglFIKTQZKKaXQZKCUUgpNBkoppdBkoJRSCk0GSiml0GSglFIKTQZKKaXQZKCUUgr4/wHQ3HDYhNeivAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "stacked = generative_model(GG, ll, rr)\n",
    "plt.plot(coverage_range, stacked, \"y-\")\n",
    "plt.plot(coverage_range, P, \"c-\")\n",
    "plt.gca().set_xlim(0, 200)\n",
    "plt.gca().set_ylim(0, 6e8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 327,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17 1.6 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
      "Iteration 0\n",
      "5.494310964329976 1.0919973352285304 [1.06101992 1.12914282 1.00002572 1.26510821 1.08946602 1.00000112\n",
      " 1.11158054 1.11541641 1.08148784 1.06800824 1.00382618 1.18306524\n",
      " 1.32743102 1.00000063 1.10276557 0.35976497] 8.277614013282319e+18\n",
      "Iteration 1\n",
      "5.665578024204201 1.094375731966334 [1.17812868e+00 1.21136296e+00 1.29554695e+00 1.31110330e+00\n",
      " 1.42241367e+00 1.16358248e+00 1.00001316e+00 1.43430056e+00\n",
      " 1.49517939e+00 1.00000304e+00 1.40629178e+00 1.56614123e+00\n",
      " 1.00056621e+00 1.18411783e+00 1.00000025e+00 1.88786903e-06] 8.277613957158651e+18\n",
      "Iteration 2\n",
      "5.564181489867448 1.0000000504874837 [1.24557659e+00 1.00000341e+00 1.34124576e+00 1.00623726e+00\n",
      " 1.26467132e+00 1.00046423e+00 1.16769721e+00 1.70042858e+00\n",
      " 1.57527458e+00 1.34831008e+00 1.60387361e+00 1.92244218e+00\n",
      " 1.20687501e+00 1.13188391e+00 1.16557496e+00 2.97162373e-09] 8.277613905638575e+18\n",
      "Iteration 3\n",
      "5.52620664056175 1.073717675549752 [8.31689248e-01 1.26720700e+00 1.00273525e+00 1.32398473e+00\n",
      " 1.00081529e+00 1.39569124e+00 1.00003281e+00 2.43378132e+00\n",
      " 1.08921640e+00 1.03316192e+00 2.93957665e+00 2.63781366e+00\n",
      " 1.89110750e+00 1.00000297e+00 1.12060717e+00 7.99938161e-15] 8.277613807010061e+18\n",
      "Iteration 4\n",
      "5.472712233523419 1.1499625974225056 [1.31357456e+00 1.13715500e+00 1.18919759e+00 1.00000568e+00\n",
      " 1.16813275e+00 1.06727921e+00 1.20172817e+00 1.00000759e+00\n",
      " 1.00000079e+00 1.12016833e+00 5.46716441e+00 4.17443488e+00\n",
      " 1.40371321e+00 1.13137632e+00 1.01614974e+00 1.36073487e-19] 8.277613697296394e+18\n",
      "Iteration 5\n",
      "5.725557362419597 1.1731469930475615 [8.70675864e-01 1.03738069e+00 1.13881347e+00 1.14321737e+00\n",
      " 1.01971259e+00 1.00009718e+00 1.00000654e+00 1.00974535e+00\n",
      " 1.21873944e+00 1.12682862e+00 1.05558232e+01 3.37841763e+00\n",
      " 1.00006682e+00 1.00068246e+00 1.26984215e+00 4.48665174e-20] 8.277613530328574e+18\n",
      "Iteration 6\n",
      "5.666723551871968 1.0445313277902137 [8.21510059e-01 1.10003985e+00 1.00029713e+00 1.00117686e+00\n",
      " 1.00018435e+00 1.03288757e+00 1.00002137e+00 1.01480067e+00\n",
      " 1.04682811e+00 1.21662447e+00 2.10356192e+01 1.03016417e+00\n",
      " 1.18999984e+00 1.01978669e+00 1.02127979e+00 4.22820902e-20] 8.27761317875098e+18\n",
      "Iteration 7\n",
      "5.909219102262106 1.3480906388967921 [6.90020457e-07 4.16504430e-01 1.23025520e+00 1.98004561e+00\n",
      " 3.92315776e-01 5.76752094e-01 1.39023290e+00 1.00002594e+00\n",
      " 1.00064885e+00 1.51930503e+00 8.43321381e+01 1.00000004e+00\n",
      " 6.08082038e-06 8.50995648e-01 1.00205926e+00 1.60767975e-20] 8.277610653088488e+18\n",
      "Iteration 8\n",
      "6.143688335749655 1.0062059563365535 [5.32251590e-07 7.51264427e-01 4.22188026e+00 8.20980048e-03\n",
      " 1.39456878e-01 1.26118496e-04 2.96280475e-01 1.15224077e+00\n",
      " 4.68545835e-01 1.11086454e+00 3.31079967e+02 9.00628734e-01\n",
      " 4.34820301e-06 1.09262523e-02 1.00113253e+00 2.52755620e-21] 8.277600350264582e+18\n",
      "Iteration 9\n",
      "5.888988460225725 1.6666197293472829 [1.50249115e-06 1.17224692e-01 6.78489105e+00 2.66646796e-03\n",
      " 6.87814111e-03 1.63434576e-06 1.42013242e-01 4.25159207e-03\n",
      " 4.34261253e-01 1.07835515e+00 1.58006539e+03 7.10943082e-08\n",
      " 4.87668206e-10 1.16530292e-02 2.46572705e+00 3.70414340e-21] 8.277550038611214e+18\n"
     ]
    }
   ],
   "source": [
    "# Optimize all variables at once\n",
    "def optimize_func(lambda_, rho, G):\n",
    "    def f(arg):\n",
    "        g = np.array(arg[:-2], dtype=np.float64)\n",
    "        l = arg[-2]\n",
    "        r = arg[-1]\n",
    "        if any(x < 0 for x in g):\n",
    "            return np.inf\n",
    "        if l <= 1:\n",
    "            return np.inf\n",
    "        if r <= 1:\n",
    "            return np.inf\n",
    "        return func(l, r, g)\n",
    "    res = minimize(f, tuple(list(G) + [lambda_, rho]), method='nelder-mead')\n",
    "    score = func(res.x[-2], res.x[-1], res.x[:-2])\n",
    "    return res.x[-2], res.x[-1], res.x[:-2], score\n",
    "\n",
    "G0 = np.ones(16)\n",
    "l0 = 17\n",
    "r0 = 1.6\n",
    "\n",
    "print(l0, r0, G0)\n",
    "def run_optimization(round):\n",
    "    ll, rr, GG = l0, r0, G0\n",
    "    for i in range(round):\n",
    "        print(\"Iteration\", i);\n",
    "        ll, rr, GG, score = optimize_func(ll, rr, GG)\n",
    "        print(ll, rr, GG, score)\n",
    "        \n",
    "run_optimization(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sequencing error and PCR duplicates (not used)\n",
    "\n",
    "Assume we have a kmer that has a small error rate $\\epsilon$ to generate an erroneous kmer when duplicating. We simulate 8 rounds of PCR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import choice, random\n",
    "from jcvi.formats.fasta import rc\n",
    "\n",
    "def generate_kmer(K=21):\n",
    "    s = ''.join([choice('ATCG') for x in range(K)])\n",
    "    return min(s, rc(s))\n",
    "\n",
    "def mutate_kmer(kmer, mu=0.01):\n",
    "    s = ''.join([(choice('ATCG') if random() < mu else k) for k in kmer])\n",
    "    return min(s, rc(s))\n",
    "\n",
    "def evolve(kmer, mu=0.01, start_coverage=1, end_coverage=30):\n",
    "    pool = [kmer] * max(start_coverage, 1)\n",
    "    while len(pool) < end_coverage:\n",
    "        pool += [mutate_kmer(kmer, mu=mu) for kmer in pool]\n",
    "    return pool[:end_coverage]\n",
    "        \n",
    "def summarize(pool):\n",
    "    from collections import Counter\n",
    "    for k, v in Counter(pool).items():\n",
    "        print(k, v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CTGCCCTGCTTTTATACAACA\n",
      "CTGCCCTGCTTTTATACAACA\n",
      "CTGCCCTGCTTTTATACAACA\n",
      "CTGCCCTGCTTTTATACAACA\n",
      "CTGCCCTGCTTTTATACAACA\n",
      "CTGCCCTGCTTTTATACAACA => CTGCCCTGCTTTTATCCAACA\n",
      "CTGCCCTGCTTTTATACAACA\n",
      "CTGCCCTGCTTTTATACAACA\n",
      "CTGCCCTGCTTTTATACAACA\n",
      "CTGCCCTGCTTTTATACAACA\n"
     ]
    }
   ],
   "source": [
    "starting_kmer = generate_kmer()\n",
    "for i in range(10):\n",
    "    dup_kmer = mutate_kmer(starting_kmer)\n",
    "    if dup_kmer == starting_kmer:\n",
    "        print(starting_kmer)\n",
    "    else:\n",
    "        print(starting_kmer, \"=>\", dup_kmer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CTGCCCTGCTTTTATACAACA 27\n",
      "CTCCCCTGCTTTTATACAACA 2\n",
      "CTACCCTGCTTTTATACAACA 1\n"
     ]
    }
   ],
   "source": [
    "pool = evolve(starting_kmer)\n",
    "summarize(pool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a mutation model which takes two parameters: \n",
    "\n",
    "- $\\mu$ probability of per-base error when duplicating\n",
    "\n",
    "- $r$ PCR duplication round, this is done with `start_coverage` and `end_coverage`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import lru_cache\n",
    "from collections import Counter\n",
    "import sys\n",
    "\n",
    "@lru_cache(maxsize=None)\n",
    "def sequencing_model(K=21, mu1000x=10, rho100X=20, target_coverage=30):\n",
    "    N = 10000\n",
    "    diff = Counter()\n",
    "    same = Counter()\n",
    "    total = 0\n",
    "    for _ in range(N):\n",
    "        if i % 100 == 0: print(\"Sample\", i, file=sys.stderr)\n",
    "        starting_kmer = generate_kmer(K=K)\n",
    "        pool = evolve(starting_kmer, mu=mu1000x/1000, \n",
    "                      start_coverage=target_coverage * rho100X // 100, \n",
    "                      end_coverage=target_coverage)\n",
    "        for pkmer, v in Counter(pool).items():\n",
    "            if pkmer == starting_kmer:\n",
    "                same[v] += 1\n",
    "            else:\n",
    "                diff[v] += 1\n",
    "            total += v\n",
    "    return same, diff, total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbd11e4f810>]"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEGCAYAAACHGfl5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAXr0lEQVR4nO3df5BdZ33f8fcXabVaR/iuSBNWdrFZQ4hI0thOXLlOHJTgUByjSUkwrRuQfzROdsqkNc4Uk1BI0nFmQIIEe4IpS+1gLCgYO8ENAofMBGgHVbYHg2VMqgDygomlSwhoVyyspJX07R/n3OFqfbX3rPbHuat9v2Z2zu5zznPuV3dW97PPc35FZiJJWtmeVXcBkqT6GQaSJMNAkmQYSJIwDCRJwOq6CzgdEXGMIsgO1V2LJC0jZwMnMvMZn/2xHE8tjYgTQDQajbpLkaRlY2JiAiAz8xmzQstyZAAcajQajfHx8brrkKRlY3BwkImJiY4zKh4zkCQZBpIkw0CShGEgScIwkCRhGEjSstKcbLL7G7tpTjYXdL/L9dRSSVpxduzZwcjOEfpW9TF9fJrRLaNsvXDrguzbkYEkLQPNySYjO0eYOjbFoSOHmDo2xcjOkQUbIRgGkrQMjB0co29V30lta1atYezg2ILs3zCQpGVgeP0w08enT2o7evwow+uHF2T/hoEkLQND64YY3TLKwOoBGv0NBlYPMLpllKF1Qwuy/+V6o7px700kaSVqTjYZOzjG8PrhOQdBeW+iicwcnLnOs4kkaRkZWje0YKOBdk4TSZIMA0lSxTCIiP6I2BYR+yNiKiIeiogrKvT7o4jIDl8Le+mcJGleqh4zuBt4FXAb8FXgeuDBiNicmbsr9B8Bvt/289QcapQkLbKuYRARm4BrgJsz87ay7R7gCWAb8JIKr/ORzPTUH0nqUVWmia4GpoE7Ww2ZeRi4C7g8IjZU2EdExNkREadXpiRpMVUJg4uBvZk5OaP9ESCAiyrs4ylgApiIiD+PiOfMrUxJ0mKqcsxgA/B0h/YD5fKcWfoeBP4MeAg4CryU4vjBz0TEpZl5pFOniOg2pdTosl6SNAdVwmAA6PShfbhtfUeZefuMpvsj4gngDuBa4H9UKVKStLiqTBNNAf0d2te2rZ+L91CcWXTKU1Mzc3C2L4opJ0nSAqkSBgcopopmarXtn8sLZuYJimknjxtIUo+oEgaPARsjYt2M9kvL5Z65vGBE9AHPA741l36SpMVTJQzuB/qAG1sNEdEP3ADsysz9Zdt5EbGxvWNE/EiH/b2BYorpk6dbtCRpYXU9gJyZD0fEfcD28pqCfcB1wPkUVyK33ANspjjdtOXrEfFhigvUjgC/RHEl82eB/7kQ/wBJ0vxVvR3FtcCt5XI98DhwVWbu6tLvg8DPA68G1gBfK/fz1sw8djoFS5IWng+3kaQVYraH23gLa0mSYSBJMgwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CSRMUwiIj+iNgWEfsjYioiHoqIK+b6YhHxiYjIiLht7qVKkhZL1ZHB3cDNwAeAm4ATwIMRcVnVF4qIVwAvmWuBkqTF1zUMImITcA1wS2bekpnvBV4KPAVsq/IiEbEGeCewfR61SpIWSZWRwdXANHBnqyEzDwN3AZdHxIYK+7gJGADecTpFSpIWV5UwuBjYm5mTM9ofAQK4aLbOETEEvAV4U2Z+/7SqlCQtqtUVttkAPN2h/UC5PKdL/7cCf09xvKGSiBjvskmj6r4kSd1VCYMB4EiH9sNt6zsqjzdcC2zOzJx7eZKkpVAlDKaA/g7ta9vWP0NEBHA78BeZ+dm5FJWZg7OtL0cOjg4kaYFUCYMDFFNFM7Xa9p+i368Bm4A3RcTzZ6w7u2z7ZmZ2DBNJ0tKpcgD5MWBjRKyb0X5pudxzin7nlfv/FDDW9gVwQ/n95jlVK0laFFVGBvcD/wW4EbgNiiuSKT7Qd2Xm/rLtPOCszNxb9vsY8LUO+/sosJPi1NTPz6d4SdLC6BoGmflwRNwHbC+vKdgHXAecD1zftuk9FH/pR9lvX7ntSYpDCezLzAfmW7wkaWFUGRlAcUbQreVyPfA4cFVm7lqswiRJSyeW4xmfETHeaDQa4+PdLkeQJLUMDg4yMTEx0emMTW9hLUkyDCRJhoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSRV0pxssvsbu2lONusuZVGsrrsASep1O/bsYGTnCH2r+pg+Ps3ollG2Xri17rIWlCMDSZpFc7LJyM4Rpo5NcejIIaaOTTGyc+SMGyEYBpI0i7GDY/St6jupbc2qNYwdHKuposVhGEjSLIbXDzN9fPqktqPHjzK8frimihZHpTCIiP6I2BYR+yNiKiIeiogrKvR7TUR8KiKaEXEkIr4WEe+LiPPnX7okLb6hdUOMbhllYPUAjf4GA6sHGN0yytC6obpLW1CRmd03ivgQ8CrgNuCrwPXAJcDmzNw9S7/twAZgD/Ad4Hzgt4FVwE9n5mlNukXEeKPRaIyPj59Od0mas+Zkk7GDYwyvH162QTA4OMjExMREZg7OXNc1DCJiE/AwcHNm3la2rQWeAPZn5kvmUkxE/AzwKPCGzHzHXPq27cMwkKQ5mi0MqkwTXQ1MA3e2GjLzMHAXcHlEbJhjPV9v1TXHfpKkRVIlDC4G9mbm5Iz2R4AALuq2g4h4TkT8aERcAryvbP7bOVUqSVo0VS462wA83aH9QLk8p8I+vgz8cPn9t4HfycxPn2rjiOg2/9Oo8JqSpIqqhMEAcKRD++G29d38OvBDwEbgtcCzK1UnSVoSVcJgCujv0L62bf2sMvP/lN8+GBEPAE9ExGRmvusU2896PKEcOTg6kKQFUuWYwQGKqaKZWm375/KCmTlGcTbRa+bST5K0eKqEwWPAxohYN6P90nK55zRedwD/speknlElDO4H+oAbWw0R0Q/cAOzKzP1l23kRsbG9Y0T8yMydRcTPUpyB9Og86pYkLaCuxwwy8+GIuA/YXl5TsA+4juJq4uvbNr0H2ExxumnL1yPiI8AXgUngJ4H/AHwXuHUh/gGSpPmr+jyDayk+vK8F1gOPA1dl5q4u/e4Afhl4JXAWxfGHjwC3lscOJEk9oNK9iXqNt6OQpLmb7+0oJElnOMNAkmQYSJIMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgSaJiGEREf0Rsi4j9ETEVEQ9FxBUV+v16RNwbEWMR8f2I2BsRb4+IxvxLlyQtlKojg7uBm4EPADcBJ4AHI+KyLv3eC7wY2AH8Z+CT5XJXRKw9nYIlSQtvdbcNImITcA1wc2beVrbdAzwBbANeMkv3qzPzMzP29yjw/nKfd59W1ZKkBVVlZHA1MA3c2WrIzMPAXcDlEbHhVB1nBkHpo+XyxdXLlCQtpiphcDGwNzMnZ7Q/AgRw0Rxfc6hc/tMc+0mSFknXaSJgA/B0h/YD5fKcOb7mG4HjwF+eaoOIGO+yDw9AS9ICqjIyGACOdGg/3La+koj4DeA3ge2Zua9qP0nS4qoyMpgC+ju0r21b31VE/ALFcYaPA2+ZbdvMHOyyr3EcHUjSgqkyMjhAMVU0U6ttf7cdRMSFwF8BjwP/LjOPV65QkrToqoTBY8DGiFg3o/3Scrlnts4R8QLgr4F/BF6Rmd+bc5WSpEVVJQzuB/qAG1sNEdEP3ADsysz9Zdt5EbGxvWNEDAF/Q3GR2ssz0zOIJKkHdT1mkJkPR8R9wPbymoJ9wHXA+cD1bZveA2ymON205a+BC4DtFNckXN62bl9m7p5f+ZKkhVDlADLAtcCt5XI9xdz/VZm5q0u/C8vlLR3WvR8wDCSpB0Rm1l3DnEXEeKPRaIyPd7scQZLUMjg4yMTExESnMza9hbUkyTCQJBkGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgSaJiGEREf0Rsi4j9ETEVEQ9FxBUV+m2KiHdHxKMRcTQicv4lS5IWWtWRwd3AzcAHgJuAE8CDEXFZl35XAb9Vfr/vdAqUJC2+rmEQEZuAa4BbMvOWzHwv8FLgKWBbl+7/HTg7M38W+OR8i5UkLY4qI4OrgWngzlZDZh4G7gIuj4gNp+qYmd/MzKl5VylJWlRVwuBiYG9mTs5ofwQI4KIFr0qStKRWV9hmA/B0h/YD5fKchSunEBHjXTZpLPRrStJKVmVkMAAc6dB+uG29JGkZqzIymAL6O7SvbVu/oDJzcLb15cjB0YEkLZAqI4MDFFNFM7Xa9i9cOZKkOlQJg8eAjRGxbkb7peVyz8KWJElaalXC4H6gD7ix1RAR/cANwK7M3F+2nRcRGxelSknSoup6zCAzH46I+4Dt5TUF+4DrgPOB69s2vQfYTHG6KQARcT6wtfxxU9n25vLnPZn5sfn+AyRJ81flADLAtcCt5XI98DhwVWbu6tJvuOzXrvXz+wHDQNKsmpNNxg6OMbx+mKF1Q3WXc8aKzOV377iIGG80Go3x8W6XI0haznbs2cHIzhH6VvUxfXya0S2jbL1wa/eO6mhwcJCJiYmJTmdsegtrST2pOdlkZOcIU8emOHTkEFPHphjZOUJzsll3aWckw0BSTxo7OEbfqr6T2tasWsPYwbGaKjqzGQaSetLw+mGmj0+f1Hb0+FGG1w/XVNGZzTCQ1JOG1g0xumWUgdUDNPobDKweYHTLqAeRF8nKO4DcbMLYGAwPw5C/VFKv82yiheMB5JYdO+CCC+DKK4vljh11VySpi6F1Q1z2vMsMgkW2ckYGzWYRAFNt99UbGIAnn3SEIGlFcGQAxdRQ38lnJrBmTdEuSSvcygmD4WGYPvnMBI4eLdolaYVbOWEwNASjo8XUUKNRLEdHnSKSJFbSMYMWzyaStELNdsyg6o3qzhxDQ4aAJM2wcqaJJEmnZBhIkgwDSZJhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQNIpNCeb7P7GbpqTzbpL0RKoFAYR0R8R2yJif0RMRcRDEXFFxb7nRsRHImI8Ig5FxAMR4X2jpR62Y88OLrj9Aq784JVccPsF7NjjUwHPdJXuWhoRHwJeBdwGfBW4HrgE2JyZu2fptw74PPBs4E+BY8DNQAIXZebB0yp6PnctlTSr5mSTC26/gKljP3gq4MDqAZ686UkfPbnMzetJZxGxCbgGuCUzb8nM9wIvBZ4CtnXp/jrghcBVmfn2zHwn8K+BcylCYeVqNmH37mKpntAr0yJ11zF2cIy+VSc/FXDNqjWMHfSpgGeyKtNEVwPTwJ2thsw8DNwFXB4RG7r0fSgzv9DWdy/wt8C/Pa2KzwQ7dhTPY77yymK5wyF43XplWqQX6hheP8z08ZOfCnj0+FGG1zu7eyarEgYXA3szc3JG+yNAABd16hQRzwJ+Gvhch9WPAC+KiLNO0Xd8ti+gUaHu3tRswsgITE3BoUPFcmTEEUKNmpNNRnaOMHVsikNHDjF1bIqRnSNL/pd5r9QxtG6I0S2jDKweoNHfYGD1AKNbRp0iOsNVebjNBuDpDu0HyuU5p+j3HKC/bbuZfaPc974KNZw5xsagr68IgZY1a4p2H7pTi9a0SPsceWtaZCk/AHulDoCtF27lZS94GWMHxxheP2wQrABVwmAAONKh/XDb+lP143T6djq40W5Zjw6Gh2H65CE4R48W7apFr0yL9EodLUPrhgyBFaTKNNEUxV/4M61tW3+qfpxm3zPX0BCMjsLAADQaxXJ01FFBjXplWqRX6tDKVGVkcIBiOmemVtv+U/T7DsWo4FR9k85TSGe+rVvhZS8rpoaGhw2CHtAr0yK9UodWniph8BhwU0Ssm3EQ+dJyuadTp8w8ERFfpLgeYaZLga9k5vfnVO2ZZGjIEOgxvTIt0it1aGWpMk10P9AH3NhqiIh+4AZgV2buL9vOi4iNHfr+q4i4uK3vj1Ncp3DfPGuXJC2QriODzHw4Iu4DtpfXFOwDrgPOp7gSueUeYDPFWUIt7wZ+C/hERPwJxRXIv0sxPfTOhfgHSJLmr8o0EcC1wK3lcj3wOMVVxbtm65SZ342IX6T44H8LxUjk08DrM/Pbp1u0JGlhVbo3Ua/x3kSSNHez3ZtouYbBCSAajeV5qYEk1WFiYgIgM/MZx4uXaxgco5hyOlR3LfPQSrKJWqvoHb4fP+B7cTLfj5PN5/04GziRmc84RLAsw+BMUF5F3fVq65XC9+MHfC9O5vtxssV6P3zSmSTJMJAkGQaSJAwDSRKGgSQJw0CShGEgScLrDCRJODKQJGEYSJIwDCRJGAaSJAyDJRUR/zIi7oiIv4uI70XEUxHx4Yh4Yd219YKIuCUiMiIeq7uWupS/Ix+PiIMRMRkReyLi+rrrqkNE/FhE3BsR/1D+f/m7iPi98rG7Z6yI2BARb4uIT0fEd8v/E794im1/NSI+HxGHy8+TP4yIqg8tO8lpddJpeyPw8xTPf34cGAJ+B/hCRGzKzP9XZ3F1iogh4M3A9+qupS4R8SvA/wI+Q/FkwGngRcDzaiyrFhFxLvAIxW2a3wV8B/gF4K3ATwJb66tu0f04xWfFVyk+J36u00bl78sDwKeA/wT8C+APgH9W/jwnnlq6hCLi54DPZebRtrYfA74IfDgzr6+rtrpFxN3AeRSj1cHMvKjeipZWRDSAL1P8HtxUdz11i4g3Am8Dfiozv9TWfj/wb4CzMnO6rvoWU0Q8G1iTmd+OiFcCHwV+KTM/M2O7LwGHgU2Zebxs+2Pg94GNmfmVubyu00RLKDP/b3sQlG1fAb4EvLiequoXEZuA1wK/W3ctNfoNYJDiLzsi4tkREfWWVKuzy+U3Z7Q3KUZMx5e2nKWTmd/t9oz4iPgJ4CeA0VYQlN5N8bn+qrm+rmFQs/I//HOBf6q7ljqU//4/A96fmSv2WAHwy8Be4KqI+AbFU/y+U84dr6q3tFr873J5V0RcGBHPi4jXANcD2zLzRH2l9YSLy+Xn2hszcz/wD23rK/OYQf1eA5wL/Ne6C6nJtRR/4byy7kJq9kKKYwN3A9uBLwBbKOaO1wKvr62yGmTm30TEW4A3Ab/atuoPMvPWmsrqJRvK5YEO6w4A58x1h4ZBjSJiI3AH8FlgR83lLLlybvRtwNsys9Mv9UqyDlgP/F5mbivb/jIi1gGvi4g/zsyVNnocoziY/lHg28ArgP8WEd/KzPfUWVgPGCiXRzqsOwycNdcdGgY1Kc+e+ThwEHj1Ch32vhk4Cvxp3YX0gKly+aEZ7R8EXg1sAj6xpBXVKCKuAUaBF5VTH1CE47OAd0TEvZl5sL4Ka9f6fel0mu3atvWVecygBuWZIw8CDeDlmdmsuaQlFxEbKKY+7gCeGxHPj4jnU/wiryl/Xl9jiUutNTKaecC09fNKei8AXgc82hYELX8F/BBw4dKX1FNavy8bOqzbAMx837oyDJZYRKwFPkZx/viWzPz7mkuqy3OBNcA2iumA1telFGdWjVHMl68Uj5bLc2e0//Ny+a0lrKUXPBfodOC8r1yu9FmN1skWl7Q3RsQ5FL8zcz4ZwzBYQuVZIfcCl1FMDT1Uc0l1GgN+rcPXl4Cvld/fU1dxNbivXP5mq6E80+pGigvxVtrvypeBSyLiBTPa/z3FaaWPL31JvaO89mIv8Nszzjb7j8AJ4C/mus+Vnq5L7U8ozoz4GPCciHht27rJzHygnrKWXmZOUFw9eZKIeD1wbCW9FwCZ+WhE3AP8fkT8KPB5igOmLwduycxDtRa49N4O/AqwKyJaVyBvKdvek5n/WGdxiy0i3lx+27r+aGtEXA6MZ+a7yrY3UEybfTIi7gV+iuKOBqOZ+eU5v6ZXIC+diPgMsPkUq7+emc9fump6U/kerbgrkAEiYg3FbSiuo7hVyZPAOzNztNbCalJejPhHFOfM/zDFaPJ9wNtnXGh1xomIU30wn/Q5UV6h/IcUofEt4M+BWzPz2Jxf0zCQJHnMQJJkGEiSDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJOD/A6AT4lwUUfk7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# rho100X is the diversity of starting material, 1 = 1% of the target coverage\n",
    "mm_same, mm_diff, mm_total = sequencing_model(mu1000x=8, rho100X=50, target_coverage=10)\n",
    "\n",
    "x, y = zip(*mm_same.items())\n",
    "z = [_x * _y / mm_total for _x, _y in zip(x, y)]\n",
    "plt.plot(x, z, \".\", color='g')\n",
    "\n",
    "x, y = zip(*mm_diff.items())\n",
    "z = [_x * _y / mm_total for _x, _y in zip(x, y)]\n",
    "plt.plot(x, z, \".\", color='r')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
